Computer model calibration based on image warping metrics: an application for sea ice deformation

Arctic sea ice plays an important role in the global climate. Sea ice models governed by physical equations have been used to simulate the state of the ice including characteristics such as ice thickness, concentration, and motion. More recent models also attempt to capture features such as fractures or leads in the ice. These simulated features can be partially misaligned or misshapen when compared to observational data, whether due to numerical approximation or incomplete physics. In order to make realistic forecasts and improve understanding of the underlying processes, it is necessary to calibrate the numerical model to field data. Traditional calibration methods based on generalized least-square metrics are flawed for linear features such as sea ice cracks. We develop a statistical emulation and calibration framework that accounts for feature misalignment and misshapenness, which involves optimally aligning model output with observed features using cutting edge image registration techniques. This work can also have application to other physical models which produce coherent structures.

Sea ice is frozen sea water that insulates the warmer ocean from the colder atmosphere in winter and provides a partial barrier to heat, moisture, and momentum transfer between the ocean and atmosphere, making it an important component in the Earth's energy balance. The ice cover waxes and wanes seasonally and at its maximum extent covers around 7% of the Earth's surface and close to 12% of the Global Ocean Weeks [2010]. Cracks in the ice cover, called leads, occupy 1-2% of the ice cover in winter but account for half of the heat flux from the ocean to the atmosphere Badgley [1961].
A constitutive model, the elastic-decohesive model [Schreyer et al., 2006] has been developed to describe crack formation in sea ice. The important aspect of this model is that the existence as well as the orientation and the opening or closing of a crack is directly represented.
Parameters defining the material properties for the elastic-decohesive model were often chosen based on experiments [Schulson, 2004] and some in situ observations [Coon et al., 1998]. There has never been a more detailed study of the appropriateness of these parameters and whether other choices would better fit observations.
With new high resolution satellite data now available, unprecedented opportunities to improve our understanding of the interaction between dynamics and thermodynamics of sea ice, and to appraise new models are now possible. With this imagery, we are able to resolve long linear features in the pack ice that evolve with time that can generally be associated with fractures or leads. These leads are derived from the small-scale, time-varying deformation of the ice cover observed using the RADARSAT Geophysical Processing System (RGPS) [Kwok, 2001].
Thus, a fresh assessment, through calibration of the model parameters, is appropriate. The focus of this paper is to develop metrics suitable for comparing simulation and observations of these linear features in sea ice, and to use these metrics to calibrate the elastic decohesive model. Calibration is the traditional statistical method to infer unknown input parameters for a model of a physical process from the corresponding observed output. The Bayesian framework O'Hagan, 2001, Higdon et al., 2004] is a natural choice for calibration where the inference is based on the posterior distribution of the model parameters given the observed outputs. Here, the likelihood term in the posterior will include the non-linear physical forward model, via, the elastic-decohesive model, which is computationally expensive to solve numerically. Moreover, the likelihood needs to be evaluated at each iteration of the Markov Chain Monte Carlo (MCMC) sampling scheme, thus making the Bayesian calibration method computationally very challenging. Alternatively, a statistical emulator can be constructed using the simulation runs of the elastic-decohesive model at some suitably chosen design points.
The emulator can then act as a surrogate, providing fast predictions of the simulation output at unsampled input values, with corresponding measures of uncertainty O'Hagan, 2002, Conti et al., 2009]. In practice, a Bayesian modularization approach [Bayarri et al., 2007, Liu et al., 2009] is often used, in which a Gaussian process (GP) based emulator is fit only once, then the fitted emulator replaces the physical model inside the likelihood for Bayesian calibration. We adapt this in our proposed method.
There are two major challenges in applying the traditional emulation-calibration framework to linear features such as sea ice cracks. (1) In the traditional framework, the output from the elastic-decohesive model would be modeled as a realization from a GP; however, emulating the response surface with linear features using a GP is inadequate. Although the GP model is attractive for emulation because of its flexibility to fit a large class of response surfaces, it is not suitable for emulating a surface with sharp gradients and discontinuities. For model outputs with linear features, a simple GP emulator would tend to smooth out those features. (2) Generalized least-square metrics are implicitly used in the likelihood functions when fitting a GP emulator and carrying out calibration; however, such metrics are also flawed when used in the presence of sharp gradients and discontinuities. Because the linear features occupy only a small, usually lower dimensional, region of the domain, models that completely fail to predict these important features may actually score better by these metrics than those that simulate the dynamically important features albeit at a somewhat different location or orientation than observations indicate [Mass et al., 2002]. We propose a novel emulation-calibration framework, which utilizes cutting-edge image registration technique, to circumvent the above challenges in dealing with linear sea ice features. It involves optimally aligning high-dimensional model outputs to the observed features, then performing emulation and calibration based on the metrics provided from image registration.
Image registration, finding a correspondence between pixels across images, is a standard problem in image processing (see [Goshtasby, 2012] for a review). Recently, the concept of warping images for registration using metrics has become increasingly popular given the properties they provide (Beg et al. [2005], Tagare et al. [2009], Vercauteren et al. [2009], Xie and Srivastava [2016]). Image registration is usually posed as a variational problem where the objective function between two images f 1 and f 2 is of the form where | · | is the Euclidean norm, γ is the deformation, R is a regularization to ensure smooth deformations, and Γ is the space of relevant deformations. Recently, this metric has been applied to model emulation and calibration [Kleiber et al., 2014]. However, the 'metric' used in Kleiber's approach is not a proper metric in the image space under warping and can fail to properly calibrate when observation and simulation images are not properly aligned.
The underlying problem with the above formulation is that it does not provide a proper metric that is inverse consistent. Here we propose a new method using the concepts of phase and amplitude of an image. The power of this method is that we have two metrics, first proposed by Xie and Srivastava [2016], which are proper distances in the space of images. One, the phase energy, gives a measure of how much warping is required to register geometric features. The other, the amplitude energy, gives a measure of the pixel intensity differences, if any, after an optimal warping has been applied. The combination of the measures allows one to determine how well a model is reproducing geometric features or coherent structures. A low amplitude energy, for instance, implies there is a smooth deformation from one image to the other, meaning a geometric feature is likely captured but perhaps misaligned. The magnitude of the phase energy gives a measure of how misaligned the feature is. In contrast, a high amplitude energy implies the geometric features of the two images are inconsistent and are not captured by the model. Using these measures in our calibration scheme allows the emulator to better distinguish how different a model prediction truly is from an observation.
We propose to use these two informative metrics to train our emulator for use in calibration. This allows us to use fewer parameter settings for expensive model runs while providing a unique measure of model performance. The advantage of our approach is two-fold: (1) the alignment avoids penalizing predicted features that are misaligned and misshapen in calibration, and (2) emulation-calibration for scalar metrics is straightforward so that standard software can be applied directly. The rest of this paper is organized as follows, Section 2 introduces the RGPS deformation data and the sea ice model that describes the physics of sea ice deformation. Section 3 introduces image warping and the emulation-calibration framework based on the metrics provided by the warp. We illustrate the usefulness of our approach via a toy example and a perfect model experiment based on the sea ice model in Section 4, and we apply the proposed method to the RGPS data and present the results in Section 5. Some final remarks on the proposed approach are included in Section 6.

Data
The fundamental quantity produced by RGPS is trajectories of marked sea-ice points. At an initial time, a set of points forming a square grid is located in a synthetic aperture radar (SAR) image. In images resulting from subsequent satellite passes, the original points are found again using area-based and feature-based tracking. The time separation between repeat observations is variable and is based on available coverage. The time interval between successive images is called the timestep. This procedure provides the trajectory of each point as these points move with the ice cover. Since the same set of points is tracked over a ice season, RGPS provides a densely sampled Lagrangian picture of the motion, similar to what would be obtained by drifting buoys. Secondary procedures in RGPS derive estimates of ice deformation. If the set of points in the original configuration is viewed as vertices of grid cells then the motion of the points determines the deformation of those cells. With this interpretation, grid quantities can be approximated that help provide a picture of the rate and type of ice deformation.
The RGPS deformation product consists of cell area changes and spatial derivatives computed using the ice displacements at the cell vertices. Each cell has a unique identifier, the time of its creation, the map coordinates of the cell center, the displacement of the cell center between this and the last observation, the timestep and the derivatives of displacement increments. Since leads are formed by cracks in the ice, the displacement associated with a lead is discontinuous. The article by Coon et al. [2007] introduced the idea of separating the displacement into the sum of a continuous and discontinuous part. A kinematic analysis of RGPS deformations was performed assuming that the predominant term was the jump discontinuity in displacement. This procedure fit a jump in displacement, or a crack, in each RGPS cell that could best account for the observed deformation. This fitting procedure provides the orientation and amount of normal and shear opening of the crack. Details on this fitting processes can be found in . Figure 5(b) shows an example of applying the procedure to a region of the Beaufort Sea to plot the magnitude of the predicted jump in displacement.
Observations and simulations of a region of the Beaufort Sea between Banks Island and Point Barrow during the period from Feb. 23 -Mar 11, 2004, are analyzed in this paper. During this period, daily observations are available on a 10 km square grid. Specifically, the domain represents a crack with a normal to the crack given by n = (cos α(s), sin α(s)) and a tangent to the crack given by t = (− sin α(s), cos α(s)) in the plane of the ice. Thus, the orientation of the crack is determined by α, the angle that the normal makes to the s 1 -axis. The size of the crack is obtained from the normal and tangential components of the jump vector, u n and u t , respectively. We threshold the jump in displacement and drop fractures smaller than 0.4 km.
Simulations using the material-point method (MPM) for this problem are described in Sulsky and  and output the same jump vector, [u], as that which can be derived from the RGPS data. To summarize, the simulations are performed on the same 10 km grid as the observations, initially with four material points per cell. Land is modeled with rigid, immobile material points. Boundaries of the rectangular region intersecting ice are given a prescribed velocity in agreement with the RGPS observations. At the start of the calculation the ice is at rest with a uniform average thickness. The thickness distribution producing the average thickness is assumed to be parabolic with a maximum at h = 3 m and equaling zero at h = 0 m and h = 6 m. Six hour wind fields from the National Center for Environmental Protection and National Center for Atmospheric Research (NCEP/NCAR) reanalysis [Kalnay et al., 1996] are used to determine the wind drag. Additionally, six hour air temperature, specific humidity, long wave flux, and short wave flux from the NCEP/NCAR reanalysis are used for the ice column temperature calculation. Details of the thermodynamic model are also in Sulsky and Peterson [2011]. The ocean currents are updated daily using the output from an ocean model (MITgcm, Marshall et al. [1997]) run independently from the ice simulation. Nominal values of sea surface temperature and salinity of −1.8 • C and 32 ppt are used throughout the domain for the calculation.
The elastic-decohesive constitutive model used in the simulations is based on two elastic parameters, Young's modulus, E, and Poisson's ratio (ν = 0.36), and five independent decohesive parameters. There are three strength parameters: tensile failure strength (τ nf ), shear failure strength (τ sf ), and compressive failure strength (f c ); plus the shear magnification factor set to 4, and the opening parameter, u 0 , set to 0.4 km. For the calibration, τ nf is chosen uniformly in the range 10-50 kPa and the shear strength in the range 10-150 kPa. These choices allow the shear strength to be both larger and smaller than the tensile strength.The compressive strength is set to 5τ nf and the Young's modulus to E = 10f c . Keeping these ratios means that u 0 would not need to be changed to maintain stability of the algorithm and avoid snapback. Note that these strength ranges are smaller than what has been seen in experiments (eg. Timco and Weeks [2010]), but are larger than what is typical for large scale simulations. The fact that sea ice on a large scale might be effectively weaker could be attributed to features such thermal bending fractures that appear on a scale of about 1 km, or other impurities and weaknesses that would not show up on a laboratory scale, but are smaller than typical simulation scales.
Fractures are initialized in the simulations using the kinematic analysis of the RGPS data over one day (23 February). Depending on deformation in a cell the calculated jump in displacement can reflect an opening mode (u n > 0) or a closing mode (u n < 0). An opening mode is initialized as a crack with the amount of opening prescribed by the kinematic analysis.
The discontinuities in a closing mode are initialized as cracks with the calculated direction and tangential jump in displacement, but with zero normal jump in displacement. After the initialization, the simulation is run for 15 more days. The initial fractures may open or close, and other fractures may be created and evolve. The resulting predicted values of u n , u t and α are compared to the kinematic analysis of RGPS data.

Statistical Method
We propose a novel emulation and calibration framework based on proper distance metrics defined on the image space, this allows us to properly account for the geometric differences between model outputs and observations. In this section, we first introduce the concept of phase and amplitude distances provided by image warping, where the first one quantifies the amount of geometric deformation required to correct for misaligned and misshapen model errors, while the latter quantifies the remaining model error after optimal alignment. We then describe how to construct an emulator and carry out calibration based on these two scalar metrics. This simplifies the high-dimensional image model outputs to scalar responses and significantly reduces the complexity of the statistical emulator and parameter inference procedure.

Image Warping
Our approach to image warping is to use the concept of phase and amplitude of an image that was introduced by [Xie and Srivastava, 2016]. This work has been motivated by the the same concept developed in the work on shape analysis of objects and phase-amplitude modes in functional data (see Srivastava and Klassen [2016], Marron et al. [2015], and Tucker et al. [2013] and references therein for a review of the concept). A deformation is defined to be a change of only the phase of an image (where the pixels are located) and not the amplitude. This view point then sets up the registration problem.
Let the set of all images be F = {f : D → R n |f ∈ C ∞ (D)} of the images in the domain D that map to pixel values in R n . Here we require n ≥ dim(D) for the theory presented in this section. An example of a suitable image is that of an RGB image. In this case, D = R 2 and n = 3, the space of red, green, and black intensities at a given s 1 , s 2 ∈ D. In the case of a gray scale image, where f (s 1 , s 2 ) : D ⊂ R 2 → R, one may instead consider an analogous image using the gradient of f such that F (s 1 , s 2 ) = (∂ s 1 f, ∂ s 2 f ) or a combination of intensity and gradient such as F (s 1 , s 2 ) = (f, ∂ s 1 f, ∂ s 2 f ). We discuss our particular choices for image type in Section 4.2. Also let Γ = Diff + (D) be the set of of all orientation-preserving diffeomorphisms that preserve the boundary of D. It can be shown that Γ is a group that contains an identity and inverse element. The registration of two images f 1 and f 2 is expressed as finding γ * such that (1) is not a proper Riemannian metric; specifically, it lacks the symmetry property.
To overcome this problem Xie and Srivastava [2016] introduced the q-map, where for an f ∈ F, we have a mapping Q : F → L 2 for any s ∈ D, where a(s) = ||Jf (s)|| V and Jf (s) is the Jacobian matrix of f at s. The term a(s) is known as the generalized area multiplication factor of f at s. The use of the q-map provides a distance that is a proper metric on the space of images for the registration problem For proofs of the metric properties of (3) the reader is referred to Xie and Srivastava [2016]. This registration problem can then be solved using gradient descent.
The gradient descent is computed iteratively over Γ, where at each iteration the directional derivative is taken on the tangent space of Γ at the identity element (T γ id (Γ)). The directional derivative is computed using an orthonormal basis expansion on T γ id (Γ) and details are given in Xie and Srivastava [2016]. The current implementation of the gradient decent scheme is sensitive to step size and the number of basis functions used to represent γ * in the orthornomal expansion. These can be tuned, however, for the specific problem.
This registration problem then gives rise to the definition of amplitude distance between two images We can then also define a distance on the amount of deformation required to register the images.
Since the space Γ does not have a simple Riemannian geometry or one cannot be found through a transformation as in [Tucker et al., 2013] we will use an extrinsic approach. We will use the q-map of the warping function and the resulting L 2 metric. This gives us a nice metric for computing the distance, but further statistical analysis becomes hard, such as principal component analysis, due to the fact that the inverse mapping is not available in closed form.
We will call this the phase distance and it is defined as where q γ = Q(γ). This distance measures how far the computed warping is from the identity warping. We will use these two distances in the calibration solution discussed later on.

Statistical Emulator
The sea ice model represents the formation and evolution of leads(cracks) resulting from the complex interactions among input processes. Due to the computational cost, we have a limited number of model runs. A standard approach is to replace the computer model by a fast surrogate statistical model often called an emulator. Typically, the first step is to obtain model outputs at a set of design points, then construct a statistical model such as a Gaussian process to interpolate these model outputs [Oakley and O'Hagan, 2002, Higdon et al., 2004, Conti et al., 2009.
Instead of interpolating model outputs, we use an emulator to interpolate the amplitude and phase distances, d a and d p , provided by image warping. We let Y(θ) ∈ F denote the model output corresponding to input parameter θ ∈ Θ ⊂ R d for an integer d ≥ 1, then the dis- We construct a Gaussian process emulator η(θ, d) with mean µ = (µ a , µ p ) T and covariance function C(·, ·) for the distance metrics d(θ). For the mean we use a linear function µ k = x T β k for both k = a and p, where x are the covariates and defined as the input parameter θ plus an intercept term, and β k is a d + 1-dimensional vector of regression coefficients. The covariance function models the dependence between every pair of metrics in terms of their distance in the parameter space. For example, we have used the exponential covariance function in our simulated example and data application. It is a commonly-used, separable, and stationary covariance function with the following form, where ρ k1 , . . . , ρ kd are the range parameters, σ 2 k is the partial sill, and τ 2 k is the nugget. The range parameters control the rate of the correlation decay with distance for each dimension of the input space. Here we use ξ k = (ρ k , σ 2 k , τ 2 k ) to denote the parameters for the statistical emulator, to distinguish them from the physical model parameters θ. We consider independent Gaussian processes for amplitude and phase distance since these two metrics are uncorrelated and hence we set Cov (d a (θ), d p (θ )) = 0. Other covariance funtions can be used here; we suggest performing cross-validation to select the appropriate covariance function.
The emulator is then constructed based on the metrics evaluated at the chosen design points using Latin hypercube sampling (LHS) [McKay et al., 1979]. Let X be an N × (d + 1) covariate matrix and d k = (d k (θ 1 ), . . . , d k (θ N )) T , k = a, p be the vectors of amplitude and phase distance metrics evaluated at the design points. The Gaussian process emulator described above specifies d a and d p as the following independent multivariate normal random variables: We estimate the emulator parameters β k and ξ k by maximizing the likelihood given by the above model and denote the resulting maximum likelihood estimators (MLE's) asβ k andξ k . For any new input parameter setting θ ∈ Θ, the emulator predicts its corresponding distance metrics through the conditional normal distribution of d a (θ) | d a and d p (θ) | d p , using the MLE'ŝ β k andξ k as plug-in estimates. This predictive process serves as the surrogate model in the following section on model calibration.

Model Calibration
We use Z d = (d a , d p ) T = (0, 0) T as the observation for calibration, since we use the observed image Z as the reference, and registering an image to itself results in zero amplitude distance and zero phase distance. We assume that the data model for Z d is where θ * is the 'best' parameter setting that yields d(θ * ) = 0. The data-model discrepancy = ( a , p ) T is assumed to follow a truncated normal distribution with the following density: with the discrepancy standard deviation ψ k > 0 for k = a, p and the probability and cumulative density functions φ and Φ for the standard normal distribution. We have chosen a truncated normal distribution with lower bound zero since we assume that the distance metrics for model outputs cannot take lower values than those for the observational data Z.
We choose to adapt the standard Bayesian calibration framework for parameter inference as it provides more detailed information on input parameter uncertainty, in the form of posterior densities [Kennedy and O'Hagan, 2001]. Each input parameter in θ receives a uniform prior with lower and upper bounds determined by its physical limits. Each discrepancy variance parameter ψ 2 k receives an informative inverse gamma prior with a shape parameter a k = 20 and a scale parameter b k = a k d 2 k,10 th for k = a, p, where d k,10 th denotes the 10 th quantile of the distance metrics d k obtained at the design points. This prior encourages the resulting full conditional posterior for the standard deviation ψ k to have its mode around d k,10th , so that roughly 10 percent of the existing model runs at the design points are covered by the 95% credible interval. Informative priors are chosen for the variance parameter because we have only one observation for calibration; this inverse problem is non-identifiable without strong priors. The posterior samples for θ, ψ 2 a and ψ 2 p are obtained from Markov chain Monte Carlo using the Metropolis-Hastings algorithm.

Experiments
We conduct two experiments to illustrate the usefulness of the proposed method. In the first experiment, we generate model outputs from applying a diffeomorphism on a template that represents a feature similar to a crack. The second experiment is based on the sea ice model outputs, where one run of the model is randomly chosen to be the synthetic observation. The experiments establish proof of concept, the first in a simplified setting and the second with a complex geometry that is consistent with the real data application. The experiments also establish the applicability of this method to discrete data sets, even though the theory itself is developed in a space of smooth functions. In both experiments, we compare calibration results using the Euclidean distance, phase distance, amplitude distance, and both phase and amplitude distance combined.

Diffeomorphism Warp
We illustrate the calibration method using a toy example with geometric features similar to the ice opening features resulting from ice deformation. The mathematical model used here applies a boundary preserving differmorphism to a 2-dimensional surface on the unit square, with input parameters θ = (θ 1 , θ 2 ) to define the map. This allows us to evaluate our calibration method on a data set where simulation differs from observation by an actual diffeomorphism, an ideal case. Let s = (s 1 , s 2 ) denote the coordinate of an image pixel. The diffeomorphism W : R 2 → R 2 maps a 2-dimensional coordinate to a new coordinate that is also 2-dimensional.

(9)
We generated 50 model runs at a set of design points chosen using Latin hypercube sampling. Figure 1a shows the template image for generating the model runs. For each design point, the template is deformed to create a warped image based on equation (9). As an illustration, Figure 1 shows the model outputs for four different input parameter settings. The branch located at the bottom right corner of the template image is stretched to the right when θ 1 is positive and compressed to the left when negative, while a similar deformation action is applied to the vertical direction controlled by θ 2 . We also generate a synthetic observation (Figure 1b) from the diffeomorphism warp using input parameter (0.3, 0.1). Therefore, all the model runs here have misaligned and misshapen features compared to the synthetic observation.
We compute the deformation metrics for each model run using the synthetic observation as the reference. In this case our image is of the form F (s 1 , s 2 s 2 ) is the pixel intensity. We calculate the gradient using finite difference approximations.
The optimal deformations for the four example runs and the warped model runs are shown in Figure 1. We see that all four warped model runs are indistinguishable to the synthetic observation visually. This indicates the deformation has been applied correctly and estimated well.
We fit independent GP emulators to each metric then perform Bayesian calibration as described in Section 3.3. In our MCMC implementation, we sample three 30,000 posterior samples at different initial values with 15,000 burn-in; MCMC convergence is assessed using the Gelman-Rubin diagnostic [Plummer et al., 2006]. Calibration results using different metrics are presented in Figure 2. It is notable that the posterior distributions from calibration using the amplitude or Euclidean alone are very flat, providing little information on the input parameter, while the posterior using the phase distance peaks near the true input parameter value, and the calibration using the combined metrics is similar to using phase dis- The results confirm that the proposed method is able to properly account for misaligned and misshapen model outputs. Because the model runs and the synthetic observation for this toy example are all warped from the same template, all model runs have the same amplitude distance except for small numerical differences due to interpolation in generating the images.
This means there exists a γ * that exactly maps each image to the synthetic observation, therefore, the phase distance quantifies the amount of warping in order to correct for misaligned and misshapen features. The calibration using the phase distance alone suggests the input parameter setting that produces the correct geometric features which are least misaligned and misshapen. This is particularly useful for model development and evaluation as a measure of a model's ability to produce correct coherent structures, a geometric feature that persists in a flow field. Moreover, for models which produce highly localized coherent structures, hurricanes for instance, the Euclidean metric may completely disqualify a model run which produced the correct hurricane characteristics just in the wrong location. In this case the model itself may be well calibrated with location error coming from external sources, such as error in an initial condition taken from observations.

Sea Ice Model Experiment
In the previous toy example our model runs were generated by an actual diffeomorphism, however in most geophysical applications model output comes from a physical numerical model. These model outputs corresponding to different parameters may not be an exact diffeomorphism warp of each other. Here we perform a perfect model experiment for the sea ice model described in Section 2. In this perfect model experiment we take the observation to be the model output for a chosen set of parameters. This experiment demonstrates, through the recovery of the chosen model parameterization using the calibration scheme, that the method works well for physical models whose output is not generated from a diffeomorphism.
We produce 80 models runs for a set of To carry out the analysis, we compute the optimal warping function for each model run.
Assessing the metrics corresponding to the model runs, we notice that the phase distances are left skewed. As a result, we compare cross-validation performance on GP emulators fitted to the original and the log transformed distances; results indicate that data transformation is necessary for this experiment to ensure the surrogate model has good prediction performance, for instance, small mean square prediction error and reasonable coverage.
Bayesian calibration and the necessary convergence diagnostic are then carried out; the posterior density functions are shown in Figure 3. The marginal posterior mode and credible interval for different metrics are summarized in Table 1

Application to Sea Ice Model and Sattelite Data
The proposed method is applied to the RGPS data discussed in Section 2. We warp each model output to RGPS data using the same choice for step size and the number of basis functions to represent γ as in the sea ice experiment in Section 4.2. MCMC convergence has been assessed using the Gelman -Rubin diagnostic [Plummer et al., 2006]  and discarding the first 15,000 for burn-in. Calibration results are presented in Figure 4. The posterior modes for the tensile and shear strength parameters from using the combined metrics are 25.45 (11.601, 48.063) and 75.71 (60.514, 95.42), respectively, with the credible intervals given in parenthesis. To establish that the posterior mode represents a viable parameterization, we run the elastic-decohesive model and assess the output with input parameters equal to the posterior mode, shown in Figure 5. In contrast to the perfect model case above, we cannot know whether there is a set of parameters which can exactly reproduce the RGPS data, in fact likely not. To quantify the performance of the model with the calibrated parameters we compare both the Euclidean and q-map distances (i.e. L 2 distance between q-maps of the images) between the model runs and the RGPS data, for both the calibrated parameter settings and the original parameter settings chosen from physical considerations in . We also apply our warping algorithm using each model output to measure the amplitude and phase distances between the runs and the RGPS data. Before preforming any warping we find that the parameter settings used in  give q-map and Euclidean distances of 23.9 and 46.3, respectively, between the model run and observations compared to that of 18.2 and 42.4 for the model run using the calibrated parameters. The amplitude distances for both model runs are similar at 11.2 for the original parameters and 11.9 for the calibrated parameters.
The most striking difference comes from the phase distances at 5.4 for the original parameters and 1.1 for the calibrated parameters. This suggests much less warping was required to align geometric features when comparing the calibrated model run to the observations. We believe the calibration based on both amplitude and phase provide a more reliable inference on input parameters than using Euclidean alone based on the perfect model study where the truth is known. In the instance shown in Figure 3 the Euclidean distance is less informative with a very diffuse posterior and a mode which misses the true value. Further, the center mode in Figure 4 encompasses that chosen from physical considerations in Sulsky . In the Appendix we show the case where our image was taken to be f (s) = (|[u]|(s), α(s)).

Discussion
The calibration method presented is the first of its kind based on the image warping metrics.  is of importance. Further, the metric itself is a proper distance provided the optimal γ is found.
Currently, we are using gradient descent, local minima are a possibility. The gradient decent scheme can be very sensitive to step size and the number of basis functions chosen to represent γ. In this work, we manually tuned the step size and number of basis functions to produce an optimal result. Future work would include looking at other optimization methods. The proposed method also reduces the high-dimensional model outputs into scalar responses, such that We illustrated the usefulness of the proposed method through a toy example and a perfect model experiment using the elastic-decohesive model. We used 50 design points for the first experiment and 80 for the latter. While the number of design points required depends on the complexity of the response surface and the number of active inputs influencing the response, the sample size we used here was sufficient for our experiments. A rule of thumb for selecting the number of design points is 10d, where d is the dimension of the input space [Loeppky et al., 2009, Harari et al., 2018. While the results in our experiments show that the proposed method performs well, we note that we have not added observation error when generating the synthetic observation. This can be a major caveat where a more rigorous investigation should be conducted. It remain future work to study how observation error could affect warping and calibration results.
The ultimate application of this method should be to real world observations. The use of image warping techniques is motivated by a need to assess sea ice model output consisting of fractures. For climate modeling, it is not necessary to get these features exactly right. In fact, it is not likely that we would be able to predict the cracks exactly since the forcing that drives the motion in the model is too uncertain. Thus, we desire a metric that can compare data and model favorably if the predicted cracks can be warped to the observed cracks that is, if the predicted cracks are somewhat misaligned or misshapen. In our study using the satellitederived RGPS data of sea ice deformation we are able to obtain informative posteriors using the combination of amplitude and phase metrics, and the posterior mode seems to represents a viable parameterization, when compared to that chosen from physical considerations in Sulsky   Magenta color "+" indicates the input parameters corresponding to the synthetic observation.
Bottom: RGPS data application: Posterior density function of input parameters using both amplitude and phase, amplitude only, and phase only.


extfieldsofstudy:
Mathematics
Computer Science

text:
An Exact No Free Lunch Theorem for Community Detection

A precondition for a No Free Lunch theorem is evaluation with a loss function which does not assume a priori superiority of some outputs over others. A previous result for community detection by Peel et al. (2017) relies on a mismatch between the loss function and the problem domain. The loss function computes an expectation over only a subset of the universe of possible outputs; thus, it is only asymptotically appropriate with respect to the problem size. By using the correct random model for the problem domain, we provide a stronger, exact No Free Lunch theorem for community detection. The claim generalizes to other set-partitioning tasks including core/periphery separation, $k$-clustering, and graph partitioning. Finally, we review the literature of proposed evaluation functions and identify functions which (perhaps with slight modifications) are compatible with an exact No Free Lunch theorem.

INTRODUCTION
A myriad of tasks in machine learning and network science involve discovering structure in data. Especially as we process graphs with millions of nodes, analysis of individual nodes is untenable, while global properties of the graph ignore local details. It becomes critical to find an intermediate level of complexity, whether it be communities, cores and peripheries, or other structures. Points in metric space and nodes of graphs can be clustered, and hubs identified, using algorithms from network science. A longstanding theoretical question in machine learning has been whether an "ultimate" clustering algorithm is a possibility or merely a fool's errand.
Largely, the question was addressed by Wolpert (1996) as a No Free Lunch theorem, a claim about the limitations of algorithms with respect to their problem domain. When an appropriate function is chosen to quantify the error (or loss), no algorithm can be superior to any other: an improvement across one subset of the problem domain is balanced by diminished performance on another subset. This is jarring at first. Are we not striving to find the best algorithms for our tasks? Yes-but by making specific assumptions about the subset of problems we expect to encounter, we can be comfortable tailoring our algorithms to those problems and sacrificing performance on remote cases.
As an example, the k-means algorithm for k-clustering is widely used for its simplicity and strength, but it assumes spherical clusters, equal variance in those clusters, and similar cluster sizes (equivalent to a homoscedastic Gaussian prior). Figure 1 shows the degraded performance on problems where these assumptions are violated.
To prove a No Free Lunch theorem for a particular task demands an appropriate loss function. A No Free Lunch theorem was argued for community detection (Peel et al., 2017), using the adjusted mutual information function (Vinh et al., 2009). 1 However, the theorem is inexact. A No Free Lunch theorem relies on a loss function which imparts generalizer-independence (formally defined below): one which does not assume a priori that some prediction is superior to another. The loss function used in the proof is only asymptotically independent in the size of the input. We present a correction: by substituting an appropriate loss function, we are able to claim an exact version of the No Free Lunch theorem for community detection. The result generalizes to other set-partitioning tasks when evaluated with this loss function, including clustering, k-clustering, and graph partitioning.

COMMUNITY DETECTION
A number of tasks on graphs seek a partition of the graph's nodes that maximizes a score function. Situated between the microscopic node-level and the macroscopic graph-level, these partitions form a mesoscopic structure-be it a core-periphery separation, a graph coloring, or our focus: community detection (CD). Community detection has been historically ill-defined (Radicchi et al., 2004;Yang et al., 2016), though the intuition is to collect nodes with high interconnectivity (or edge density) into communities with low edge density between them. The task is analogous to clustering, in which points near one another in a metric space are grouped together.
To assess whether the formulation of community detection matches one's needs, one performs extrinsic evaluation against a known ground truth clustering. This ground truth can come from domain knowledge of real-world graphs or can be planted into a graph as a synthetic benchmark. After running community detection on the graph, some similarity or error measure between the computed community structure and the correct one can be computed.
No bijection between true structure and graph Unfortunately, ground truth communities do not imply a single graph-and vice versa. Peel et al. (2017) go as far as to claim, "Searching for the ground truth partition without knowing the exact generative mechanism is an impossible task." We can imagine the following steps for how problem instances are created, given that we have N = |V | nodes: 1. Sample (true) partition T ∈ Ω; 2. Generate graph G from T by adding edges according to the edge-generating process g.
where Ω is our universe: the space of all partitions of N = |V | objects. Given a graph G = (V, E), we can imagine multiple truths T i ∈ Ω that could define its edge set E by different generative processes g i : Ω → Γ, where Γ is the set of all graphs with N nodes. Peel et al. (2017) give a proof that extends from this simple example: Imagine that T 1 partitions the N nodes into N components (the N-partition), and T 2 partitions them into 1 component (the 1-partition). Let g 1 exactly specify the number of edges between each pair of communities, such that g 1 (T 1 ) is G with probability 1. Similarly, let g 2 be an Erdős-Rényi model such that g 2 (T 2 ) is G with nonzero probability. (Peel et al. (2017) note that this is easily extended to graphs with more nodes.) We thus have two different ways to create a single graph; how can a method discern the correct one, without knowledge of g?
Community detection is then an ill-posed inverse problem: Use a function f : Γ → Ω to produce a clustering C = f (G), which is hopefully representative of T (Peel et al., 2017, Appendix C). 2 The function f is not a bijection, so there isn't a unique T represented in the given graph. Our algorithm f must encode our prior beliefs about the generative process g to select from among candidates. For this reason, we must hope that the benchmark graphs that we use are representative of the generative process for graphs in our real-world applications. That is, we hope that our benchmark domain matches our practical domain.
Other set-partitioning tasks While the remainder of this work focuses on community detection, our claims are relevant to other set-partitioning tasks. Notable examples are clustering (the vector space analogue to community detection), graph k-partitioning, and k-clustering. Metadata about the nodes and edges, such as vector coordinates, are used to guide the identification of such structure, but the tasks are all fundamentally set-partitioning problems. They can also have different universes Ω-the latter tasks have a smaller universe than does community detection, for a given graph G: They consider only partitions with a fixed number of clusters.

NO FREE LUNCH THEOREMS
The No Free Lunch theorem in machine learning is a claim about the universal (in)effectiveness of learning algorithms. Every algorithm performs equally well when averaging over all possible input-output pairs. Formally, for any learning method f , the error (or loss) L of the method f , summed over all possible problems (g, T ) equals a loss-specific constant Λ(L): defining the edge-generative process g and partition T as above. This loss is thus generalizerindependent. To reduce loss on a particular set of problems means sacrificing performance on others-"there is no free lunch" (Wolpert, 1996;Schumacher et al., 2001). Judiciously choosing which set to improve involves making assumptions about the distribution of the data: as we've mentioned, k-means is a method for k-clustering which works well on data with spherical covariance, similar cluster sizes, and roughly equal class sizes. When these assumptions are violated, performance suffers and overall balance is achieved.

COMMUNITY DETECTION AS SUPERVISED LEARNING
We follow Peel et al. (2017) in framing the task of community detection (CD) as a learning problem. While recent algorithms, e.g. , have introduced learnable parameters to community detection algorithms the CD literature's algorithms are by and large untrained. These untrained algorithms encode knowledge of the problem domain in prior beliefs. We note that our work and Peel et al. (2017) straightforwardly handle both of these cases.
In general supervised machine learning problems, we seek to learn the function that maps an input space X to an output space Y. We consider problem instances as sampled from random variables over each, so our goal is to learn the conditional distribution p(Y | X). In the process of training on a dataset D, we develop a distribution over hypotheses q which are estimates of the distribution p.
In the case of most community detection algorithms, our input space is the set of graphs on N nodes Γ, and the output space is ω. There is no training data: D = ∅. All of our prior beliefs about p must be encoded in the prior distribution Pr(q). That is, the model itself must contain our beliefs about the definition of community structure. Only from the encoded Pr(q) and an observed x ∈ X (our graph G) do we form our point estimate of the true distribution p (Peel et al., 2017). However, in the case of trainable CD algorithms, we encode our beliefs in the posterior distribution Pr(q | D).

LOSS FUNCTIONS AND A PRIORI SUPERIORITY
How should we evaluate an algorithm's predictions? Classification accuracy won't cut it: When comparing to the ground truth, there are no specific labels (e.g. no notion of a specific "Cluster 2")-only unlabeled groups of like entities. We settle for a measure of similarity in the groupings, quantifying how much the computed partition tells us about the ground truth.
A popular choice of measure is the normalized mutual information (NMI; Kvalseth, 1987) between the prediction and the ground truth. While this measure has a long history in community detection, its flaws have been well-noted (Vinh et al., 2009;Peel et al., 2017;McCarthy and Matula, 2018;McCarthy et al., 2019). It imposes a "geometric" structure upon the universe Ω, 3 so something as simple as guessing the trivial all-singletons clustering outperforms methods that try at all to find a mesoscopic-level structure (McCarthy et al., 2019). The property which NMI lacks is generalizerindependence.
The property of generalizer-independence is defined by the generalization error function, an expectation of the loss E[L | p, q, D]. To satisfy this property, the generalization error must be independent of the particular true value T . This is best expressed by Equation 1.
The adjusted mutual information (AMI, defined in section 3) (Vinh et al., 2009) is a proposed replacement for NMI which does not impose a geometric structure upon the space. Unfortunately, this benefit is not fully realized when the expectation is computed over a space Ψ ⊂ Ω. For the Ψ used in Peel et al. (2017), the expected AMI across all problems is only asymptotically generalizer-independent as the graph size grows-it is within some diminishing amount of error ε(N) of generalizer-independence, as proven by Peel et al. (2017). Peel et al. (2017) frame community detection in the style of learning algorithms, letting them prove a No Free Lunch theorem for community detection. They note that the claim holds for "an appropriate choice of . . . L"-specifically a loss function L that is generalizer-independent-but their chosen loss function is not fully generalizer-independent. They also consider a stricter property than generalizerindependence: homogeneity. With a homogeneous loss function, the distribution of the error (not just its expectation) is identical, regardless of the ground truth. A measure which deviates from homogeneity may have this deviation bounded by a function of the number of vertices (the graph order). Lemma 1 (Peel et al., 2017). Adjusted mutual information (AMI) is a homogeneous loss function over the interior of the space of partitions of N objects, i.e., excluding the 1-partition and the N-partition. Including these, AMI is homogeneous within 1 B N . 4

PREVIOUS RESULT: APPROXIMATE NO FREE LUNCH THEOREM
Wolpert (1996) gives a generalized No Free Lunch theorem, which assumes a homogeneous loss. Theorem 1 (Wolpert, 1996). But this choice of measure (AMI) is not, in fact, homogeneous over the entire universe Ω (Lemma 1). A strategy that guesses either of the non-interior (i.e., boundary) partitions-the 1-partition or Npartition-will yield a higher-than-average reward. There is indeed a negligible amount of free lunch-a free morsel, if you will.  Correcting this is arguably a pedantic demand, for two reasons: 1. The fraction 1 B N converges to 0 superexponentially as N increases. 2. The deficiency is only present when T is one of the trivial partitions. Otherwise, AMI as used is exactly homogeneous. But the trivial partitions reflect a lack of any mesoscopic community structure.
Nevertheless, we'd like to see a tight claim of generalizer independence. To do this, we must select the proper random model, a sample space for a distribution.
AMI adjusts NMI by subtracting the expected value from both the numerator and the denominator, shown in blue: where I is the mutual information, maximized when the specific clustering C equals the ground truth T . By inspecting Equation 2, we see that AMI's value is 1 (the maximum) when C = T , 0 in expectation, and negative when the agreement between C and T is worse than chance.
Subtly hidden in this equation is the decision of which distribution to compute the expectation over. For decades, this distribution has been what Gates and Ahn (2017) call M perm : all partitions of the same partition shape 5 as C or T . For example, if C partitioned 7 nodes into clusters of sizes 2, 2, and 3, then we would compute the expected mutual information over all clusterings where one had cluster sizes of 2, 2, and 3. McCarthy et al. (2019) argue that M perm is inappropriate. To use this random model assumes that we can only produce outputs within that restricted space, when in actuality Ω is the set of all partitions of N nodes. Furthermore, during evaluation, we hold our ground truth fixed, rather than marginalizing over possible ground truths. Were we to instead consider a distribution over T s, we would add noise from other possible generative processes which yield the same graph from different underlying partitions. In our average, we might be including scores on ground truths that better align with our notions of, say, core-periphery partitioning. For this reason, we take a one-sided expectation-over C, holding T fixed. The one-sided distribution over all partitions of N nodes is called M 1 all Gates and Ahn (2017). This distribution is what we use for our AMI expectation, giving a measure denoted as AMI 1 all , which is recommended by McCarthy et al. (2019). It takes the form AMI 1 all (C, T ) The differences between M all and M perm are illustrated in Figure 2 under |V | = 3. We will now show that substituting M all for M perm , hence using AMI 1 all , allows for an exact No Free Lunch theorem.

AN EXACT NO FREE LUNCH THEOREM
We strengthen the No Free Lunch theorem for community detection given by Peel et al. (2017) by using an improved loss function, AMI 1 all , for community detection. Our proof does not distinguish the "boundary" partitions (the two trivial partitions) from the "interior" partitions (the remainder). It is entirely agnostic toward the particular ground truth T , which is exactly what we need. We improve the previous result by moving from M interior (which excludes the boundary partitions) to M all .

GENERALIZER-INDEPENDENCE OF AMI 1 all
Lemma 2. AMI 1 all is a generalizer-independent loss function over the entire space M all of partitions of N objects.
Proof. Like Peel et al. (2017), we must show that the sum of scores is independent of T : where Ω is the space of all partitions of N objects. Unlike Peel et al. (2017), we take the AMI expectation over all B N clusterings in Ω using the random model M 1 all (Gates and Ahn, 2017). To prove our claim about Equation 4, we note that denominator of AMI 1 all is a constant with respect to C (Equation 3), so we can factor it out of the sum and restrict our attention to the numerator. This is because the max-term in the denominator is the constant log N (Gates and Ahn, 2017) and the expectation term for a given T is independent of the particular C. Having factored this out, we will now prove Equation 4 by the stronger claim: To prove Equation 5, we separate the summation's two terms: The expectation is uniform over the universe Ω, 6 so we can apply the law of the unconscious statistician, then push the constant probability out, to get Because the inner sum is independent of any particular C , the outer sum is a sum of constantsone for each element in Ω. We can now express Equation 5 as follows, where the reciprocals straightforwardly cancel out: This equivalence implies that Equation 4 is true.
The proof is valid without loss of generality vis-à-vis the distribution-that is, as long as the AMI expectation is computed uniformly over the problem universe Ω, AMI is a generalizer-independent measure. This stipulation is relevant to tasks which assume a fixed number of clusterings-using M num -like k-clustering and graph partitioning.
Having demonstrated the generalizer-independence of AMI, we can define our loss function as, say, The loss is zero when we exactly match the true clustering and positive otherwise.
Having proven the generalizer-independence of AMI 1 all , we now turn to a more general form of the No Free Lunch theorem, which admits not just a homogeneous loss function but any generalizerindependent loss. Theorem 3 (Wolpert, 1996). For generalizer-independent loss , the uniform average over all p, There is no free lunch.") Proof. See Wolpert (1996).
Theorem 4 (No Free Lunch theorem for community detection and other set-partitioning tasks). For a set-partitioning problem with a loss function of adjusted mutual information using the appropriate random model for the task, the uniform average over all p, E [ | p, D], equals Λ( ) |Y| .
Proof. Lemma 2 proves that AMI using the appropriate random model is generalizer-independent. Applying Theorem 3 completes the proof (Peel et al., 2017).

OTHER MEASURES
AMI stemmed from a series of efforts to improve normalized mutual information (NMI). We note that six other measures, when extended to M 1 all instead of M perm , are also generalizer-independent: the adjusted Rand index (ARI; Hubert and Arabie, 1985), relative NMI (rNMI; Zhang, 2015), ratio of relative NMI (rrNMI; , Cohen's κ (Liu et al., 2018), corrected NMI (cNMI Lai and Nardini, 2016), and standardized mutual information (SMI; Romano et al., 2014). We elide the proofs because they are similar to Lemma 2. Each of the six measures satisfies the precondition for the No Free Lunch theorem when the random model matches the problem domain.
Of late, a renewed push has advocated using the adjusted Rand index (ARI; Hubert and Arabie, 1985) to evaluate community detection; in fact, ARI and AMI are specializations of the same underlying function which uses generalized information-theoretic measures (Romano et al., 2016). Every claim in the proof works for ARI, by replacing every mutual information I term with the Rand index RI.
Another line of research, focusing on improving NMI, produced rNMI Zhang (2015), rrNMI , and cNMI (Lai and Nardini, 2016). We note that rrNMI is identical to one-sided AMI when both are extended to M 1 all . Consequently, our claim above works just as well for rrNMI. Further, because we were able to ignore the denominator of AMI in our proof of Lemma 2, we can do the same for rrNMI, which gives its unnormalized variant, rNMI. This means that rNMI is a generalizer-independent measure as well, when used in the appropriate one-sided random model. The practical benefit of normalizing rNMI into rrNMI is that the normalized measure gives a more interpretable notion of success.
Additionally, Lemma 2 holds true for standardized mutual information (which is equivalent to standardized variation of information and standardized V-measure) (Romano et al., 2014), the adjusted variation of information (Vinh et al., 2009), and for Cohen's κ, advocated for CD by Liu et al. (2018). This is because each measure shares the form of AMI: an observed score minus an expectation.
Finally, to show whether cNMI is generalizer-independent under the correct random model, we must show how to specialize it into a one-sided variant, because there is room for interpretation about how this should be done, even restricting our focus to M 1 all . The expression for cNMI depends on both C and T relative to the universes that contain them. Our specialization should remove dependence on the family of T , so we arrive at the following expression after cancellation and noting that the NMI between a clustering and itself is 1: As it turns out, this quasi-adjusted measure is also generalizer-independent.
In general, we now have a recipe for generalizer-independent loss functions: They can be created by subtracting the expected score from the observed score. This recipe works whenever a uniform expectation can be well defined.

CONCLUSION
We now have a proof of the No Free Lunch theorem for community detection and clustering that is both complete and exact. We show that a corrected form of AMI, namely AMI 1 all , computes its expectation in a way that does not advantage the boundary partitions (1 cluster and N singleton clusters). Indeed, this expectation is over the entire universe of partitions Ω, rather than any proper subset, such as the historically common M perm . We affirm the claim: "Any subset of problems for which an algorithm outperforms others is balanced by another subset for which the algorithm underperforms others. Thus, there is no single community detection algorithm that is best overall" (Peel et al., 2017).
It is still possible for an algorithm to perform better on a subset of community detection problems, so we can strive toward improved results on such a subset. To hope to perform well, we must note the assumptions about the subset of problems we expect to encounter. Some work has been done on estimating network properties to select the correct algorithm for the task at hand-a coarse way of checking assumptions (Peel, 2011;Yang et al., 2016). Beyond this, though, we must clarify what the problem of community detection is; the formulation we choose will guide which subset of problem instances to prioritize and which to sacrifice.