# Is In-Domain Data Really Needed?
A Pilot Study on Cross-Domain Calibration for Network Quantization

Haichao Yu1 Linjie Yang2 Humphrey Shi1
1University of Illinois at Urbana-Champaign  2ByteDance Inc.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Post-training quantization methods use a set of calibration data to compute
quantization ranges for network parameters and activations. The calibration
data usually comes from the training dataset which could be inaccessible due
to sensitivity of the data. In this work, we want to study such a problem: can
we use out-of-domain data to calibrate the trained networks without knowledge
of the original dataset? Specifically, we go beyond the domain of natural
images to include drastically different domains such as X-ray images,
satellite images and ultrasound images. We find cross-domain calibration leads
to surprisingly stable performance of quantized models on 10 tasks in
different image domains with 13 different calibration datasets. We also find
that the performance of quantized models is correlated with the similarity of
the Gram matrices between the source and calibration domains, which can be
used as a criterion to choose calibration set for better performance. We
believe our research opens the door to borrow cross-domain knowledge for
network quantization and compression.

##  1 Introduction

With the increasing popularity of deploying neural networks on edge devices,
neural network quantization has become a widely studied topic [6, 32, 27, 43,
20, 42, 39, 15, 2, 41]. By quantizing its weights and activations to low-bit
integers, a neural network can be stored with smaller size and executed at
faster speed with less memory footprint and computational resources.

Most existing network quantization methods can be roughly divided into two
groups. One is quantization-aware training (QAT), the other is post-training
quantization (PQ). By inserting differentiable simulated quantization
operations into the network during training, QAT methods allow training losses
to back-propagate through the quantization operations. Although QAT methods
can achieve satisfactory performance in most cases, it is time-consuming and
data-hungry, which is unacceptable in some scenarios, , when the time budget
is restricted or training dataset is inaccessible. Compare with QAT, PQ has
advantages on time and data efficiency. Given a pretrained full-precision
model, PQ runs the model over a small set of calibration data to calculate
value ranges of intermediate activations. These calibration data shares the
same distribution as the training data on which the model is trained. After
this, both network activations and weights can be quantized into low-precision
integers for inference acceleration. Compared with QAT, this post-training
quantization process is highly efficient (, in a few minutes).

Figure 1: Pipeline of cross-domain calibration for post-training quantization
without using in-domain data.

However, PQ methods still require in-domain calibration data, which may be
inaccessible in some situations due to privacy or security reasons. For
example, federated learning [16] trains models using distributed private user
data that are not accessible by the developers. Without calibration data from
the same domain as the training data, the traditional PQ methods will fail.
Targeting at this problem, Cai et al.proposed a novel method termed ZeroQ that
generates synthetic data for calibration [3]. Specifically, ZeroQ utilized
BatchNorm [13] statistics from the pretrained model as supervision to optimize
randomly initialized input such that the BatchNorm statistics of the input are
similar to those from the pretrained model. ZeroQ showed its effectiveness on
various datasets including ImageNet [7] and MSCOCO [23]. Another method is to
use Generative Adversarial Networks (GAN) to create fake calibration data
[37]. The proposed framework named GDFQ is trained with BatchNorm statistics
loss, Cross Entropy classification loss and knowledge distillation loss.

One drawback of these data-free quantization methods is their high time
complexity. With common hardware setting (, on single GPU), ZeroQ requires
hundreds of back-propagation iterations to synthesize a single batch and GDFQ
needs to be trained for hundreds of epochs. On the other hand, due to the
optimization nature of these methods, hyper-parameters also need to be tuned
for satisfactory results. If in-domain training data is inaccessible and
synthetic data generation is time-consuming, a straightforward question to ask
is “Can we use real-world images from another domain for calibration?”. As
shown in Figure 1, we use calibration data from domain $B$ to calibration
model trained on domain $A$, resulting in a quantized model. The model is
further evaluated on a test set in domain $A$.

To study this problem, we carry out a large-scale benchmark with 10 different
tasks and 13 diversified image domains including nature images, low-resolution
images, ultrasound images, satellite images, etc. We find that a simple
treatment on the BatchNorm layers in the calibration procedure greatly
improves performance of cross-domain calibration, almost bridging the gap
between in-domain and out-of-domain calibration on a wide range of tasks. We
also find that the performance of quantized models is correlated with the
similarity of the Gram matrices between the source and calibration domains,
which can be used as a criterion to choose similar image domains for better
performance. Compared with synthetic data generation method ZeroQ, our
approach achieves comparable or better performance on our large-scale
benchmark. Although the study is mainly empirical, it reveals that cross-
domain data can be used for post-training calibration just as effective as in-
domain data, which could motivate the community to explore more in the
direction of model quantization and compression with cross-domain knowledge.

##  2 Related Work

###  2.1 Quantization-Aware Training

Quantization-aware training (QAT) inserts simulated quantization operation
into model forward and backward passes during training. First introduced by
Courbariaux et al. [12], binary neural networks achieved seven times faster
inference speed on GPU. Further in [32], Rastegari et al.proposed to use a
binary tensor and a float scalar to approximate the weight or activation
tensor. Parallel to binarization, multi-bit model quantization is also widely
explored. DoReFa-Net [42] was proposed to quantize both network and training
gradients. Scale-adjusted training [14] improve the model performance by
scaling activation in the network to preserve proper scale of gradients. To
bring flexibility to the model quantizer, some works also introduced learnable
quantization schemes. In [39], LQ-Net was proposed to learn the quantization
basis vectors. In [15], Jung et al.proposed a novel method to learn the
quantization intervals.

###  2.2 Post-training Quantization

Different from QAT, post-training quantization (PQ) directly computes
quantization parameters with a calibration set. In [17], Krishnamoorthi et
al.introduced PQ with layer-wise and channel-wise quantization schemes. To
improve quantization performance, Banner et al.proposed a novel per-channel PQ
scheme that analytically computes clipping thresholds and bits allocated for
each channel [2]. To mitigate the quantization error from the clipping
operation, Zhao et al.introduced a technique named OCS to split the channels
with extreme values [41].

Traditionally, PQ method requires a calibration dataset which is often sampled
from training data to estimate the clipping ranges of model activations.
However, in some scenarios, training data may be inaccessible due to data
privacy issues. In this case, quantization could be conducted without a
calibration set. Assuming features from a BatchNorm layer follow Gaussian
distribution, Nagel et al.directly used BatchNorm statistics (, running means
and variances) to estimate activation ranges [29]. The drawback of this method
is the estimated ranges may be inaccurate since the distribution of
intermediate layers may be different from Gaussian distribution. DFC and ZeroQ
proposed to synthesize calibration data by gradient descent under the
supervision of BatchNorm statistics [10, 3]. Generative Adversarial Networks
(GAN) is also utilized to generate synthetic calibration data [37], in which
the generator generates calibration data and the discriminator is the
quantized model. One problem with these optimization-based methods is the high
time complexity. Hundreds of iterations or even hundreds of epochs of back-
propagation are commonly required. In addition, much human efforts are needed
for hyper-parameter tuning such as learning rate and optimizer selection. In
contrast to these methods, we find that out-of-domain image datasets can serve
as calibration data effectively for a wide range of tasks.

###  2.3 Batch Normalization and Domain Knowledge

To adapt a deep neural network to a different domain, Li et al.proposed
Adaptive BatchNorm to update running statistics in target domain [22]. The
central assumption is that domain specific information is encoded in BatchNorm
layers. Sharing a similar spirit with this, Li et al.showed that style
transfer can be conducted by matching the BatchNorm statistics between two
images [21]. Inspired by the effect of BatchNorm statistics, we design a
BatchNorm updating method on out-of-domain calibration dataset to estimate
ranges of activations, which effectively improves the performance of cross-
domain calibration even when the training and the calibration domains have
huge appearance differences.

##  3 Cross-domain Calibration for Post-training Quantization

###  3.1 Motivation

Given a pretrained full-precision model, post-training quantization (PQ) is
widely used for inference acceleration. Most existing PQ methods require a set
of calibration data to calculate the quantization parameters for weights and
activations. Training data is normally used for calibration. However, in some
situations, the training data is not available due to privacy or security
issues. To solve this problem, Cai et al. [3] proposed to synthesize
calibration data using BatchNorm statistics as supervision. Since different
types of real images are vastly available, we would like to explore another
option: use out-of-domain real images for calibration. We conduct a large-
scale empirical study with drastically different datasets and tasks to
investigate this setting.

###  3.2 Quantization Schemes

In all our experiments, we employ layer-wise uniform quantization for both
network weights and activations [17]. To quantize a tensor $\bm{x}$ by
$k$-bit, the quantized tensor $\bm{x}_{q}$ is calculated as

$\displaystyle\bm{x}_{int}$
$\displaystyle=\mathrm{round}\left(\frac{\bm{x}}{s}\right)+z,$
$\displaystyle\bm{x}_{q}$
$\displaystyle=\mathrm{clip}\left(\bm{x}_{int},c_{l},c_{h}\right),$ (1)

where $s$ and $z$ are scale and zero point parameters, $c_{l}$ and $c_{h}$ are
lower and upper clipping thresholds. For network weight $\bm{w}$, we use a
symmetric min-max quantization scheme, which is defined as

$\displaystyle s_{w}$
$\displaystyle=\frac{\mathrm{max}(\left|w_{min}\right|,\left|w_{max}\right|)}{2^{k-1}},$
$\displaystyle z_{w}$ $\displaystyle=0,$ $\displaystyle c_{l}$
$\displaystyle=-2^{k-1},$ $\displaystyle c_{h}$ $\displaystyle=2^{k-1}-1,$ (2)

where $w_{min}$ and $w_{max}$ are minimum and maximum of $\bm{w}$. To quantize
activation $\bm{a}$, we employ an affine histogram quantization scheme, which
is also used in
PyTorch111https://pytorch.org/docs/master/torch.quantization.html. To
determine the quantization parameters, we have

$\displaystyle s_{a}$ $\displaystyle=\frac{a_{h}-a_{l}}{2^{k}-1},$
$\displaystyle z_{a}$ $\displaystyle=\frac{a_{l}(2^{k}-1)}{a_{l}-a_{h}},$
$\displaystyle c_{l}$ $\displaystyle=0,$ $\displaystyle c_{h}$
$\displaystyle=2^{k}-1,$ (3)

where $a_{l}$ and $a_{h}$ are lower and upper clipping thresholds for
$\bm{a}$. For each layer, a histogram is built to represent the activation
distribution. We search for the optimal $a_{l}$ and $a_{h}$ values such that
the quantization error is minimized with respect to $\bm{a}$. In all our
experiments, we first fold BatchNorm layers to their preceding convolutional
or linear layers before computing quantization parameters.

###  3.3 Datasets and Networks

(a) Imagenet [7]

(b) MSCOCO [23]

(c) Cityscapes [5]

(d) Pascal VOC 07 [8]

(e) Cifar100 [18]

(f) CelebA [26]

(g) CUB-200-2011 [35]

(h) Agriculture (NIR and RGB channel) [4]

(i) Kaggle Ultrasound [1]

(j) NIH Chest X-ray [36]

(k) Open Images [19]

(l) Oxford Flowers [31]

(m) SVHN [30]

Figure 2: Example images from different domains. We will use abbreviations in
the following section: IN: ImageNet, CO: MSCOCO, CS: Cityscapes, VO: Pascal
VOC 2007 CI: Cifar100, CE: CelebA, CB: CUB-200-2011, AG: Agriculture-Vision,
US: Kaggle Ultrasound, NI: NIH Chest X-ray. OI: Open Images, OF: Oxford
Flowers-101, SV: SVHN.

Our experiments span across different image domains including natural images
[19, 7, 8, 5, 23, 18], X-ray [36], ultrasound [1], and satellite images [4].
We also experiment with models in various computer vision tasks including
image classification, semantic segmentation and object detection. The datasets
and corresponding example images are illustrated in Figure 2. We train
floating-point models on all the datasets from Figure 2(a) to Figure 2(j),
where each dataset has one or more associated models. These datasets and three
additional datasets from Figure 2(k) to Figure 2(m) are used as calibration
sets. When calibration and training data have different sizes, we resize the
calibration data to the same size of the training data. One exception is when
calibrating ResNet-18 on Cifar100, calibration samples are randomly cropped.
In addition, when Agriculture-Vision data is used to calibrate models on other
datasets, only RGB channels are used. When other datasets are used to
calibrate Agriculture-Vision models, we use RGB channels to compute a
grayscale channel to be used as NIR channel.

These models to be quantized are summarized in Table 1. The evaluation metrics
are summarized as below:

  1. 1.

All classification tasks except CelebA: top-1 accuracy.

  2. 2.

CelebA: average top-1 accuracy on 40 attributes.

  3. 3.

Pascal VOC 2007: mAP.

  4. 4.

MSCOCO, Cityscapes and Agriculture-Vision: mIoU.

  5. 5.

Ultrasound: DICE.

Specifically, DICE is a widely used evaluation metric for binary segmentation,
which is defined as

$Dice=\frac{2|\bm{p}\cap\bm{g}|}{|\bm{p}|+|\bm{g}|},$ (4)

where $\bm{p}$ and $\bm{g}$ are binary prediction and ground truth tensors
respectively.

Datasets | Tasks | Models
---|---|---
Imagenet | C |  ResNet-18, ResNet-50 [11]
Inception-V3 [34]
MSCOCO | S |  FCN [28]
Cityscapes | S |  BiSeNet [38]
Pascal VOC 2007 | D |  MobileNetV2 SSD-Lite [33] [25]
Cifar100 | C | ResNet-18
CelebA | C | ResNet-50
CUB-200-2011 | C |  MMAL-Net [40]
Agriculture-Vision | S |  MSCG-Net-101 [24]
Ultrasound | S | U-net
NIH Chest X-ray | C | ResNet-50
Table 1: Training datasets and the corresponding pretrained models. C: Image
Classification. S: Semantic Segmentation. D: Object Detection.

###  3.4 Cross-domain Calibration

Figure 3: Cross-domain calibration with and without BatchNorm adjustment.
Dashed green line denotes performance of in-domain calibration. The dataset
abbreviations are defined in Figure 2. Zoom in for details.

#### A Naive Approach

First we employ a naive cross-domain calibration approach, , directly using
out-of-domain data to calibrate a model trained in another domain. The results
of 8-bit quantization are plotted as blue bars in Figure 3. The in-domain
calibration results are marked in dashed lines.

We can see that for all the datasets and models, a large number of out-of-
domain calibration results are comparable to the in-domain baselines. Take
ResNet-50 on ImageNet as an example, the gap of out-of-domain and in-domain
calibration results is very small. This observation is surprising since most
existing PQ methods assume in-domain calibration data is required to estimate
activation ranges of intermediate layers. For some models such as MobileNetV2
SSD-Lite on Pascal VOC, ResNet-50 on CelebA, U-net on Ultrasound, and
ResNet-50 on NIH Chest X-ray, there is still an accuracy gap between out-of-
domain calibration results and the baseline in-domain results. Such
discrepancy often occurs when the original and the calibration domain have a
large gap, most natural image datasets are drastically different from
ultrasound images. Motivated by this observation, we introduce an approach to
reduce the representation gap between the original and the calibration domains
for the target model, namely BatchNorm adjustment.

#### BatchNorm Adjustment

Since BatchNorm was proposed to reduce Internal Covariate Shift [13], it has
been widely applied to deep neural networks. For one neural network, different
domains will generate different BatchNorm parameters. One interesting
application is that BatchNorm parameters can be updated on a new domain for
domain adaptation [22, 21]. As stated by these work, BatchNorm layers encode
information that is specific to the domain that the model is trained on.
Inspired by this, we propose BatchNorm adjustment to adapt the models to the
out-of-domain calibration datasets. We describe our method as below:

  1. 1.

Given a model $M$ to be quantized and an out-of-domain dataset $D$.

  2. 2.

Reset running means and variances of all BatchNorm layers in $M$:
$\mu\leftarrow 0$, $\sigma\leftarrow 1$.

  3. 3.

Run $M$ on $D$ to accumulate new BatchNorm statistics.

  4. 4.

Fold BatchNorm layers to their preceding convolutional or linear layers.

  5. 5.

Calculate weight quantization parameters ($s_{w}$ and $z_{w}$ in Equation 2).

  6. 6.

Run $M$ on $D$ again to calculate activation quantization parameters ($s_{a}$
and $z_{a}$ in Equation 3).

The results are plotted in orange bars in Figure 3. For most datasets and
models, BatchNorm adjustment gets better performance than the baselines
without BatchNorm adjustment. In most cases, the results are comparable to or
even slightly better than the in-domain calibration results. On ResNet-18/50
on ImageNet, FCN on MSCOCO, ResNet-18 on Cifar100, U-net on Ultrasound and
ResNet-50 on NIH Chest X-ray, all out-of-domain calibration results are within
0.3% performance gap from the in-domain calibration results. On some tasks
such as MMAL-Net on CUB-200-2011 and MSCG-Net-101 on Agriculture-Vision, some
out-of-domain data still performs relatively lower than others. In such cases,
randomly chosen calibration datasets cannot guarantee performance similar to
in-domain calibration. In order to fix this issue, we investigate how domain
discrepancy affects calibration and how to improve the calibration results.

###  3.5 Influence of Domain Discrepancy on Calibration

Figure 4: Visualization of correlation between calibration performance and
domain discrepancy. Top-3 out-of-domain datasets with smallest domain
discrepancy from the in-domain data are marked in  red. Dashed green line
denotes performance of in-domain calibration. The dataset abbreviations are
defined in Figure 2. Zoom in for details.

In this section, we investigate the relationship between domain discrepancy
and cross-domain calibration performance. Since Maximum Mean Discrepancy (MMD)
was proposed to measure difference of sample mean in Reproducing Kernel
Hilbert Space [9], it has been widely used as a domain discrepancy measure. As
proved by Li et al. [21], matching MMD with the second order polynomial kernel
is equivalent to matching Gram matrices of feature maps. Similar to [21], we
employ mean $L_{2}$ distance between Gram matrices of feature maps to measure
the discrepancy between two image domains. Formally, domain discrepancy $D$ is
defined as

$D=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\lVert\bm{G}_{ij}^{A}-\bm{G}_{ij}^{B}\lVert_{2}^{2},$
(5)

where $\bm{G}^{A}\in\mathbb{R}^{N\times N}$ and
$\bm{G}^{B}\in\mathbb{R}^{N\times N}$ are average Gram matrices of two domains
$A$ and $B$, which are defined as

$\bm{G}^{A}=\frac{1}{|A|}\sum_{k=1}^{|A|}\bm{G}_{k}^{A}\quad\textrm{and}\quad\bm{G}^{B}=\frac{1}{|B|}\sum_{k=1}^{|B|}\bm{G}_{k}^{B},$
(6)

where $|A|$ and $|B|$ are number of samples in domain $A$ and $B$
respectively. Ignoring the domain superscript for simplicity, each element
$\bm{G}_{k,ij}$ in $\bm{G}_{k}^{A}$ or $\bm{G}_{k}^{B}$ is

$\bm{G}_{k,ij}=\sum_{m=1}^{M}\bm{F}_{k,im}\bm{F}_{k,jm},$ (7)

where $\bm{F}_{k}$ is the feature embedding of $k$-th sample, which is
normalized by mean and standard deviation over all images in this domain, and
$M$ is the embedding dimension. In our experiments, we use the output of last
conv layer in the fourth conv block conv_4_3 of VGG-16 to extract the feature
maps.

Next we study how domain discrepancy is related to calibration performance. In
Figure 4, we plot model performance vs domain discrepancy with the source
training data into scatter points. We also calculate correlation coefficient
between model performance and domain discrepancy shown in the title of each
sub-figure. First, on most of the datasets, calibration performance is
negatively correlated to gram matrix distance. The out-of-domain datasets that
have the smaller domain discrepancy with the in-domain data always achieve
comparable performance with in-domain data. We marked the top-3 out-of-domain
data with smallest domain discrepancy in red for better visualization. This
observation can be useful for cross-domain calibration, when
calibration/training data is sensitive and inaccessible but some high-level
statistics such as Gram matrix are available. One can ask the data owner to
provide a mean Gram matrix of the source dataset, which can be used to search
in the pre-built candidate pool of cross-domain datasets for the out-of-domain
calibration dataset with smallest domain discrepancy.

Calib Methods \ Datasets | IN |  CO |  CS |  VO |  CI |  CE |  CB |  AG |  US |  NI
---|---|---|---|---|---|---|---|---|---|---
R18 |  R50 |  IV3 |  |  |  |  |  |  |  |  |
FP32 |  69.76 |  76.13 |  77.46 |  60.47 |  79.10 |  0.686 |  76.50 |  91.74 |  89.63 |  54.63 |  0.730 |  0.773
In-domain |  69.39 |  75.73 |  76.88 |  60.17 |  79.09 |  0.685 |  76.30 |  91.72 |  89.36 |  53.93 |  0.738 |  0.758
ZeroQ |  69.36 |  75.56 |  76.89 |  60.17 |  78.99 |  0.681 |  76.01 |  91.59 |  87.90 |  53.34 |  0.739 |  0.718
ZeroQ-real |  69.29 |  75.57 |  76.89 |  60.18 |  78.87 |  0.681 |  75.87 |  91.62 |  89.29 |  53.29 |  0.738 |  0.675
Cross-domain |  69.47 |  75.70 |  76.93 |  60.12 |  79.00 |  0.685 |  76.36 |  91.50 |  89.19 |  53.85 |  0.733 |  0.761
Cross-domain (MS) |  69.47 |  75.70 |  76.93 |  60.08 |  79.01 |  0.684 |  76.37 |  91.50 |  89.19 |  53.38 |  0.732 |  0.761
Table 2: 8-bit quantization results with different calibration datasets. The
dataset abbreviations are defined in Figure 2. R-18, R-50 and IV3 are
Resnet-18, Resnet-50 and InceptionV3 respectively. Cross-domain (MS) uses an
average of gram matrices from multiple layers. Best results without using in-
domain data are emphasized in bold.

(a) ImageNet

(b) Open Images

(c) Ultrasound

(d) ZeroQ

Figure 5: Visualization of Gram matrices of different datasets. The feature
extraction model is VGG-16 with BatchNorm trained on ImageNet. In Column (b)
and (c), BatchNorm adjustment is applied before calculating gram matrices. In
Column (d), VGG-16 is used to synthesize the ZeroQ data. Best viewed in color.

Using the above strategy, we compare our method with in-domain calibration and
data synthesis approach ZeroQ [3] in Table 2. In addition, we also compare
with a variant of ZeroQ that we term as ZeroQ-real, short for ZeroQ with real
images. Different from ZeroQ that initializes the synthesized data with random
Gaussian values, ZeroQ-real uses real images for initialization. In all the
experiments of ZeroQ-real, ImageNet data is used for initialization. We show
performance of two variants of our method with one or multiple feature maps in
VGG-16 to compute domain discrepancy, shown as Cross-domain and Cross-domain
(MS) respectively. When multiple feature maps are used, each feature map will
generate a Gram matrix and the domain discrepancy is defined as the average of
$L_{2}$ distances of all Gram matrices. As shown in Table 2, our method
achieves comparable or even better results than both ZeroQ and ZeroQ-real.
Using one or multiple feature maps to compute domain discrepancy achieves
similar performance across different tasks. Our method using selected cross-
domain calibration data is on-par with in-domain data, proving the
effectiveness of using cross-domain calibration data for post-training
quantization.

In Figure 5, we show a visual comparison of gram matrices of different
datasets. We show Gram matrices of ImageNet, and two datasets that have the
smallest and largest domain discrepancy to ImageNet (, Open Images and
Ultrasound). In addition, we also show Gram matrix of synthetic data from
ZeroQ. We use VGG-16 with BatchNorm trained on ImageNet as feature extractor.
BatchNorm adjustment is applied when using Open Images and Ultrasound
datasets. As shown in Figure 5, Open Images shares a similar gram matrix with
ImageNet. In contrast, Ultrasound dataset has a very different Gram matrix
since it is not natural images. On the other hand, the Gram matrix of the
synthetic ZeroQ data is also drastically different from the the two natural
image datasets ImageNet and Open Images, showing less spatial correlation than
the others. Our hypothesis is that with BatchNorm statistics as the only
guidance to synthesize images is not sufficient to build spatial correlations
as in real images. The lack of spatial correlation might be a key reason that
ZeroQ does not achieve satisfactory performance on tasks such as NIH Chest
X-ray classification.

###  3.6 Visualization of Activation Ranges

Figure 6: Visualization of calculated activation ranges of different methods. Three examples are shown from left to right: ResNet-18 on ImageNet, FCN on MSCOCO and ResNet-50 on NIH Chest X-ray. In each sub-figure, upper and lower clipping thresholds are plotted separately. Calib Methods \ Datasets | IN |  CO |  CS |  VO |  CI |  CE |  CB |  AG |  US |  NI
---|---|---|---|---|---|---|---|---|---|---
R18 |  R50 |  IV3 |  |  |  |  |  |  |  |
FP32 |  69.76 |  76.13 |  77.46 |  60.47 |  79.10 |  0.686 |  76.50 |  91.74 |  89.63 |  54.63 |  0.730 |  0.773
In-domain |  62.19 |  66.80 |  65.81 |  45.97 |  52.57 |  0.612 |  74.86 |  91.00 |  87.01 |  49.11 |  0.739 |  0.711
ZeroQ |  61.63 |  66.27 |  64.76 |  52.07 |  50.21 |  0.614 |  74.18 |  90.37 |  81.98 |  48.73 |  0.741 |  0.523
ZeroQ-real |  61.80 |  66.22 |  63.87 |  52.11 |  49.78 |  0.611 |  74.43 |  90.49 |  86.54 |  47.82 |  0.742 |  0.502
Cross-domain |  62.18 |  66.82 |  65.52 |  50.90 |  53.43 |  0.608 |  74.12 |  90.02 |  87.06 |  48.37 |  0.729 |  0.712
Cross-domain (MS) |  62.18 |  66.82 |  65.52 |  50.22 |  53.87 |  0.608 |  74.44 |  90.02 |  87.06 |  48.15 |  0.729 |  0.716
Table 3: 6-bit quantization results with different calibration datasets. The
dataset abbreviations are defined in Figure 2. Best results without using in-
domain data are emphasized in bold.

In this section, we visualize the activation ranges (clipping thresholds)
calculated on different out-of-domain datasets. In Figure 6, we show three
examples: Ultrasound to calibrate ResNet-18 on ImageNet, NIH Chest X-ray to
calibrate FCN on MSCOCO, and ImageNet to calibrate ResNet-50 on NIH Chest
X-ray. In each sub-figure, we show the lower and upper clipping thresholds for
each layer. Compared with the naive cross-domain calibration method, the
proposed BatchNorm adjustment makes a closer estimation of activation ranges
to those from the in-domain data. Specifically, for ResNet-50 on NIH Chest
X-ray, BatchNorm adjustment shows significant improvement over the baseline
with original BatchNorm parameters. On the other hand, our estimation is also
better than ZeroQ showing that calibration with cross-domain data is more
robust than synthetic data in these cases.

###  3.7 Lower-bit Quantization

We also explore cross-domain calibration for 6-bit quantization. We use the
same weight and activation quantization schemes as used in 8-bit experiments.
Again, we show performance of two variants of our method with one or multiple
feature maps to compute domain discrepancy. The results are summarized in
Table 3. Best calibration results without using in-domain data are in bold. In
most cases, our proposed cross-domain method achieves comparable or better
performance than ZeroQ and its variant. Specifically, on Cityscapes and NIH
Chest X-ray, our method outperforms ZeroQ by a large margin.

##  4 Discussion

In this work, we show the feasibility of using out-of-domain data for post-
training quantization, which is different from the assumption of existing
works that in-domain calibration dataset is necessary. Data synthesis methods
are time-consuming and may not perform well in some cases as shown in our
experiments. With stable and superior performance on a wide range of tasks,
our method can be a new direction of post-training quantization when in-domain
data are not available.

There are some interesting topics to explore in future works. First, cross-
domain calibration can be improved by designing a better dataset pool
consisting of a large amount of diverse domains. Second, better domain
discrepancy measures can be explored to further improve the performance of
quantized models.

##  5 Conclusion

In this work, we explored cross-domain calibration for post-training
quantization. To study this problem, we conducted a large-scale study that
spans across various tasks, datasets and neural networks. We find that a
simple BatchNorm adjustment strategy can effectively improve the performance
of quantized models by a large margin, almost bridging the gap between cross-
domain calibration and in-domain calibration. In addition, we find that
performance with cross-domain calibration is correlated with Gram matrix
similarity between the source and the calibration domains. Therefore, Gram
matrix similarity can be used as a criterion to select calibration dataset
from a candidate pool to further improve performance. We believe our work will
motivate future research on utilizing cross-domain knowledge for network
quantization and compression.

## References

  * [1] Ultrasound nerve segmentation.  https://www.kaggle.com/c/ultrasound-nerve-segmentation/overview, 2016.
  * [2] Ron Banner, Yury Nahshan, and Daniel Soudry.  Post training 4-bit quantization of convolutional networks for rapid-deployment.  In Advances in Neural Information Processing Systems, pages 7950–7958, 2019.
  * [3] Yaohui Cai, Zhewei Yao, Zhen Dong, Amir Gholami, Michael W Mahoney, and Kurt Keutzer.  Zeroq: A novel zero shot quantization framework.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13169–13178, 2020.
  * [4] Mang Tik Chiu, Xingqian Xu, Yunchao Wei, Zilong Huang, Alexander G Schwing, Robert Brunner, Hrant Khachatrian, Hovnatan Karapetyan, Ivan Dozier, Greg Rose, et al.  Agriculture-vision: A large aerial image database for agricultural pattern analysis.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2828–2838, 2020.
  * [5] Marius Cordts, Mohamed Omran, Sebastian Ramos, Timo Rehfeld, Markus Enzweiler, Rodrigo Benenson, Uwe Franke, Stefan Roth, and Bernt Schiele.  The cityscapes dataset for semantic urban scene understanding.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3213–3223, 2016.
  * [6] Matthieu Courbariaux, Itay Hubara, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Binarized neural networks: Training deep neural networks with weights and activations constrained to+ 1 or-1.  arXiv preprint arXiv:1602.02830, 2016.
  * [7] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * [8] M. Everingham, L. Van Gool, C. K. I. Williams, J. Winn, and A. Zisserman.  The PASCAL Visual Object Classes Challenge 2007 (VOC2007) Results.  http://www.pascal-network.org/challenges/VOC/voc2007/workshop/index.html.
  * [9] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  The Journal of Machine Learning Research, 13(1):723–773, 2012.
  * [10] Matan Haroush, Itay Hubara, Elad Hoffer, and Daniel Soudry.  The knowledge within: Methods for data-free model compression.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8494–8502, 2020.
  * [11] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [12] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Binarized neural networks.  In Advances in neural information processing systems, pages 4107–4115, 2016.
  * [13] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  volume 37 of Proceedings of Machine Learning Research, pages 448–456, Lille, France, 07–09 Jul 2015. PMLR.
  * [14] Qing Jin, Linjie Yang, and Zhenyu Liao.  Towards efficient training for neural network quantization, 2019.
  * [15] Sangil Jung, Changyong Son, Seohyung Lee, Jinwoo Son, Jae-Joon Han, Youngjun Kwak, Sung Ju Hwang, and Changkyu Choi.  Learning to quantize deep networks by optimizing quantization intervals with task loss.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4350–4359, 2019.
  * [16] Jakub Konečný, Brendan McMahan, and Daniel Ramage.  Federated optimization:distributed optimization beyond the datacenter, 2015.
  * [17] Raghuraman Krishnamoorthi.  Quantizing deep convolutional networks for efficient inference: A whitepaper.  arXiv preprint arXiv:1806.08342, 2018.
  * [18] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * [19] Alina Kuznetsova, Hassan Rom, Neil Alldrin, Jasper Uijlings, Ivan Krasin, Jordi Pont-Tuset, Shahab Kamali, Stefan Popov, Matteo Malloci, Alexander Kolesnikov, Tom Duerig, and Vittorio Ferrari.  The open images dataset v4: Unified image classification, object detection, and visual relationship detection at scale.  IJCV, 2020.
  * [20] Fengfu Li, Bo Zhang, and Bin Liu.  Ternary weight networks.  arXiv preprint arXiv:1605.04711, 2016.
  * [21] Yanghao Li, Naiyan Wang, Jiaying Liu, and Xiaodi Hou.  Demystifying neural style transfer.  arXiv preprint arXiv:1701.01036, 2017.
  * [22] Yanghao Li, Naiyan Wang, Jianping Shi, Jiaying Liu, and Xiaodi Hou.  Revisiting batch normalization for practical domain adaptation.  arXiv preprint arXiv:1603.04779, 2016.
  * [23] Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick.  Microsoft coco: Common objects in context.  In European conference on computer vision, pages 740–755. Springer, 2014.
  * [24] Qinghui Liu, Michael C. Kampffmeyer, Robert Jenssen, and Arnt-Borre Salberg.  Multi-view self-constructing graph convolutional networks with adaptive class weighting loss for semantic segmentation.  In The IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, June 2020.
  * [25] Wei Liu, Dragomir Anguelov, Dumitru Erhan, Christian Szegedy, Scott Reed, Cheng-Yang Fu, and Alexander C Berg.  Ssd: Single shot multibox detector.  In European conference on computer vision, pages 21–37. Springer, 2016.
  * [26] Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.  Deep learning face attributes in the wild.  In Proceedings of International Conference on Computer Vision (ICCV), December 2015.
  * [27] Zechun Liu, Baoyuan Wu, Wenhan Luo, Xin Yang, Wei Liu, and Kwang-Ting Cheng.  Bi-real net: Enhancing the performance of 1-bit cnns with improved representational capability and advanced training algorithm.  In Proceedings of the European conference on computer vision (ECCV), pages 722–737, 2018.
  * [28] Jonathan Long, Evan Shelhamer, and Trevor Darrell.  Fully convolutional networks for semantic segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3431–3440, 2015.
  * [29] Markus Nagel, Mart van Baalen, Tijmen Blankevoort, and Max Welling.  Data-free quantization through weight equalization and bias correction.  In Proceedings of the IEEE International Conference on Computer Vision, pages 1325–1334, 2019.
  * [30] Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, and Andrew Y Ng.  Reading digits in natural images with unsupervised feature learning.  2011\.
  * [31] Maria-Elena Nilsback and Andrew Zisserman.  Automated flower classification over a large number of classes.  In 2008 Sixth Indian Conference on Computer Vision, Graphics & Image Processing, pages 722–729. IEEE, 2008.
  * [32] Mohammad Rastegari, Vicente Ordonez, Joseph Redmon, and Ali Farhadi.  Xnor-net: Imagenet classification using binary convolutional neural networks.  In European conference on computer vision, pages 525–542. Springer, 2016.
  * [33] Mark Sandler, Andrew Howard, Menglong Zhu, Andrey Zhmoginov, and Liang-Chieh Chen.  Mobilenetv2: Inverted residuals and linear bottlenecks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4510–4520, 2018.
  * [34] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2818–2826, 2016.
  * [35] C. Wah, S. Branson, P. Welinder, P. Perona, and S. Belongie.  The Caltech-UCSD Birds-200-2011 Dataset.  Technical Report CNS-TR-2011-001, California Institute of Technology, 2011\.
  * [36] Xiaosong Wang, Yifan Peng, Le Lu, Zhiyong Lu, Mohammadhadi Bagheri, and Ronald M Summers.  Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2097–2106, 2017.
  * [37] Shoukai Xu, Haokun Li, Bohan Zhuang, Jing Liu, Jiezhang Cao, Chuangrun Liang, and Mingkui Tan.  Generative low-bitwidth data free quantization.  arXiv preprint arXiv:2003.03603, 2020.
  * [38] Changqian Yu, Jingbo Wang, Chao Peng, Changxin Gao, Gang Yu, and Nong Sang.  Bisenet: Bilateral segmentation network for real-time semantic segmentation.  In European Conference on Computer Vision, pages 334–349. Springer, 2018.
  * [39] Dongqing Zhang, Jiaolong Yang, Dongqiangzi Ye, and Gang Hua.  Lq-nets: Learned quantization for highly accurate and compact deep neural networks.  In Proceedings of the European conference on computer vision (ECCV), pages 365–382, 2018.
  * [40] Fan Zhang, Meng Li, Guisheng Zhai, and Yizhao Liu.  Multi-branch and multi-scale attention learning for fine-grained visual categorization, 2020.
  * [41] Ritchie Zhao, Yuwei Hu, Jordan Dotzel, Christopher De Sa, and Zhiru Zhang.  Improving neural network quantization without retraining using outlier channel splitting.  arXiv preprint arXiv:1901.09504, 2019.
  * [42] Shuchang Zhou, Yuxin Wu, Zekun Ni, Xinyu Zhou, He Wen, and Yuheng Zou.  Dorefa-net: Training low bitwidth convolutional neural networks with low bitwidth gradients.  arXiv preprint arXiv:1606.06160, 2016.
  * [43] Chenzhuo Zhu, Song Han, Huizi Mao, and William J Dally.  Trained ternary quantization.  arXiv preprint arXiv:1612.01064, 2016.


††thanks: The authors contributed equally.††thanks: The authors contributed
equally.

# Morphological attractors in natural convective dissolution

Jinzi Mac Huang1,2<EMAIL_ADDRESS>Nicholas J. Moore3<EMAIL_ADDRESS>1\. NYU-ECNU Institute of Physics and Institute of Mathematical Sciences, New
York University Shanghai, Shanghai, 200122, China
2\. Applied Math Lab, Courant Institute, New York University, New York, NY
10012, USA
3\. Mathematics Department, United States Naval Academy, Annapolis, MD 21402,
USA

###### Abstract

Recent experiments demonstrate how a soluble body placed in a fluid
spontaneously forms a dissolution pinnacle — a slender, upward pointing shape
that resembles naturally occurring karst pinnacles found in stone forests.
This unique shape results from the interplay between interface motion and the
natural convective flows driven by the descent of relatively heavy solute.
Previous investigations suggest these structures to be associated with shock-
formation in the underlying evolution equations, with the regularizing Gibbs-
Thomson effect required for finite tip curvature. Here, we find a class of
exact solutions that act as attractors for the shape dynamics in two and three
dimensions. Intriguingly, the solutions exhibit large but finite tip curvature
without any regularization, and they agree remarkably well with experimental
measurements. The relationship between the dimensions of the initial shape and
the final state of dissolution may offer a principle for estimating the age
and environmental conditions of geological structures.

Ever-changing geological features on this planet never fail to capture our
imagination and inspire new scientific advances. Often, striking features
appear when fluid and solid interact, ranging from centimeter scale pebble
stones [1, 2] to the kilometer scale karst terrains [3, 4]. Even planetary-
scale plate tectonics are believed to have such a fluid-structure interaction
origin [5, 6, 7, 8].

Figure 1: Dissolution-induced sharpening. (a) Limestone structures form the
stone forests of Borneo (Grant Dixon). (b)-(c) Dissolution of lab-scale planar
and axisymmetric objects unveils the sharpening process; images from the same
set of experiments reported in [9]. The observed noise in the curvature
measurements results from surface impurities, like bubbles, affecting image
tracking. The final radius of curvature at the tip was measured to be 60
$\mu$m.

The direct study of geophysical structures presents unique challenges owing to
the vast range of scales, along with the limitation of only seeing the current
state. On the other hand, laboratory-scale experiments combined with judicious
physical models have proven valuable in explaining certain formations [10,
11], like the growth of icicles [12], river meandering [13, 14, 15], the
formation of stalactites and stalagmites [16, 17, 18], meteor ablation [19],
and plate tectonics [5, 6, 7, 8]. In this letter, we investigate one such
geomorphological problem, namely the formation of karst pinnacles [20, 3]. We
will demonstrate the unusual shape dynamics that result in convergence to a
morphological attractor.

Commonly seen in South Asia and the island of Madagascar [21, 22], Fig. 1(a)
shows the typical shape of the karst pinnacles that comprise stone forests.
While their origins remain unclear, studies have related such pinnacles to the
dissolution process [20, 4, 23], as many of these rocks were once immersed
under water, and the rock material is slightly water-soluble. Two questions
naturally arise: How does the rock evolve into individual pinnacles? Why does
each pinnacle exhibit the common feature of a sharp apex?

Aimed at addressing such questions, recent experiments employed lab-scale
soluble objects to recreate the stone forests purely from the perspective of
dissolution and fluid dynamics [9]. These experiments show stone forests to
manifest from a single porous, soluble block, highlighting the sharpening of
each karst pinnacle as the key to such formations. In these and other [24]
experiments, no external flow is imposed, rather the transport of relatively
heavy solute sustains a natural convective flow that drives shape evolution.

Huang et al. 2020 and Pegler & Wykes 2020 proposed a boundary-layer based
model capable of predicting sharpening [9, 25]. Notably, the model reduces to
a single integro-PDE that governs shape evolution, denoted here as the
sharpening equation (SE). Initial numerical evidence and scaling analysis of
the SE suggested shock formation and finite-time blowup of the tip curvature
[9]. Likewise, similarity solutions of a matched-asymptotic approximation
predict unbounded growth of tip curvature for certain initial conditions [25,
26].

Fig. 1 shows experimental images of dissolving planar and axisymmetric bodies
(see [9] for experimental details). Measurements of the tip curvature indeed
increase over time, as seen in Fig. 1(b) and (c), but interestingly give no
clear indication of singular behavior. To reconcile these observations,
previous studies appealed to the thermodynamic Gibbs-Thomson (GT) effect [27,
28], which regularizes the SE and limits the curvature growth. However, the
strength of the GT term used in previous simulations was at the high end of
the range estimated from physical considerations (1-10 $\mu$m) [9], thus
calling into question whether this term accurately modeled a physical effect
or was simply acting to regularize the numerics. For context, the experiments
show in Fig. 1 reach a final tip radius of 60 $\mu$m, suggesting that the GT
effect is secondary. As such, fundamental questions remain: Does the SE
support geometric shock formation? Is there a blowup in tip curvature, and, if
so, is the blowup only limited in practice by microscale thermodynamics?

Here, we resolve these and other questions by finding a class of exact
solutions to the SE in 2D and 3D that serve as attractors for the shape
dynamics. The solutions exhibit large, but finite, tip curvature, indicating
that the GT effect is not needed to regularize sharpening. Improved numerical
methods, specially tailored to the hyperbolic nature of the SE, show how
initially convergent characteristics bend to avoid crossing and eventually
straighten in pursuit of the attracting morphology. Revisited experiments
confirm the convergence to these exact solutions, thus raising the possibility
of using the solutions to infer properties of natural structures.

The model.—In accordance with Fick’s law, a soluble interface retreats with
normal velocity proportional to the gradient of the solute field
$V_{n}\propto\nabla c\cdot\mathbf{n}$ [2, 29, 30]. These dynamics can be
greatly complicated by the presence of a fluid flow, which significantly
distorts the field $c$ and alters local gradients. The flow may be forced
externally [1, 31, 32, 33, 34, 35, 36, 37] or driven by buoyancy variations
[24, 38, 9], as in the present study. The evolution of flow, solute, and body
shape are thus inextricably linked.

Due to the large Schmidt and Grashof numbers ($\text{Sc}\sim 10^{3}$ and
$\text{Gr}\sim 10^{9}$, see SI) of the pinnacle experiments, these convective
flows are confined to narrow boundary layers, enabling an explicit expression
for the 2D interface velocity [39, 9]:

$V_{n}=-a\cos^{\frac{1}{3}}\theta\,\left(\int_{0}^{s}\cos^{\frac{1}{3}}\theta\,ds^{\prime}\right)^{-\frac{1}{4}}$
(1)

where the surface tangent angle $\theta=\theta(s,t)$ is parameterized by the
arclength $s$ from the apex, as illustrated in Fig. 2(c). The constant
$a\approx 10^{-7}\mbox{ m}^{5/4}/\mbox{s}$ contains all material and fluid
properties. For simplicity, we focus on the 2D case in this Letter, with
analogous analysis for axiysymmetric (3D) objects available in the SI.

The $\theta$-$L$ formulation [40, 41, 42, 30] offers a single, scalar equation
that fully describes shape evolution:

$\frac{\partial\theta}{\partial t}=\frac{\partial V_{n}}{\partial
s}+V_{s}\frac{\partial\theta}{\partial s}.$ (2)

As above, $\theta$ represents the surface tangent angle, and the Cartesian
coordinates can easily be recovered from
$\frac{d}{ds}(x,y)=(\sin\theta,\cos\theta)$. The artificial tangential
velocity $V_{s}=\int_{0}^{s}V_{n}\partial_{s}\theta\,ds^{\prime}$ enforces an
invariant metric with respect to arclength, thereby separating $s$ and $t$ as
independent variables. Equation 2 with interface velocity Eq. 1 is the
nonlinear integro-PDE proposed in [9], here called the sharpening equation
(SE); see [25] for the Cartesian counterpart .

Previous investigations employed a finite-difference scheme to solve Eq. 2,
but with the GT regularization required to maintain numerical stability [9].
Other studies employed a matched-asymptotic expansion, but with approximation
error that may grow large with time [25, 26]. In contrast, we introduce a
method to directly propagate characteristics of Eq. 2, with no regularization
and no additional model approximation made.

To that end, consider a location $s=S^{(0)}$ on the initial geometry, with
tangent angle $\Theta^{(0)}=\theta(S^{(0)},0)$. The trajectory $S(t)$ evolves
via the ODE:

$\dot{S}(t)=\left(R\,\frac{\partial V_{n}}{\partial
s}-V_{s}\right)\Big{\rvert}_{s=S(t)},\quad S(0)=S^{(0)},$ (3)

where $R=-(\partial\theta/\partial s)^{-1}=\kappa^{-1}$ is the radius of
curvature. Combining Eqs. 2 and 3 shows that the tangent angle remains
constant along such a characteristic, $\theta(S(t),t)=\Theta^{(0)}$, thus
providing an implicit solution for any initial profile
$\Theta^{(0)}=\theta(S^{(0)},0)$. This is the essence of the method of
characteristics.

A PDE-based interpretation of Eq. 3 is also possible via implicit functions.
That is, regard $s=s(\theta,t)$, where $\theta\in(0,\pi)$ is now the
independent variable, to obtain:

$\displaystyle\frac{\partial s}{\partial t}=-\frac{\partial
V_{n}}{\partial\theta}-V_{s}\,,$ (4) $\displaystyle
V_{n}=-a\cos^{\frac{1}{3}}\theta\left(\int^{\pi/2}_{\theta}R(\theta^{\prime},t)\cos^{\frac{1}{3}}\theta^{\prime}d\theta^{\prime}\right)^{-\frac{1}{4}}\,,$
(5)

where now $R(\theta,t)=-\partial s/\partial\theta$ and
$V_{s}=\int_{\pi/2}^{\theta}V_{n}(\theta^{\prime})d\theta^{\prime}$.
Crucially, the reformulation in terms of $\theta$ implies increased numerical
tip-resolution in proportion to the sharpening. We thus solve Equations 4 and
5 numerically (see the SI for implementation details) for a class of left-
right symmetric initial conditions.

Figure 2:  Simulating dissolution-induced sharpening. (a) Evolution of the
initial shape $\theta=\operatorname{arccot}{s/\ell}$ in 2D. (b) Zooming-in
near the apex illustrates the strong sharpening effect. (c) Model schematic.
(d) Profiles of the tangent-angle $\theta(s,t)$ show a steep gradient develop
near the apex, $s=0$, consistent with (e) a tip curvature that increases by 5
orders of magnitude. (f) Characteristic curves show contours of constant
$\theta$, with the physical trajectories shown in (a) with red.

Results.—As a first numerical test, we simulate the dissolution of the initial
profile $\theta(s,0)=\operatorname{arccot}(s/\ell)$, with $\ell=1$ and $a=1$
(for other values, time could be rescaled by the factor $\ell^{5/4}/a$). As
seen in Fig. 2(a), dissolution causes the apex to sharpen as the body retreats
downwards and diminishes in size. Figure 2(b) shows a few representative
shapes at different stages of dissolution, illustrating the dramatic
sharpening effect. Figure 2(d) shows the corresponding distributions of the
tangent angle, $\theta(s,t)$. Here, a rapid change of tangent angle develops
at the tip, as is consistent with the increasing curvature
$\kappa=-\partial\theta/\partial s$ there. Indeed, the rescaled tip curvature
$\bar{\kappa}_{0}(t)=\kappa_{0}(t)/\kappa_{0}(0)$ shown in Fig. 2(e) increases
by 5 orders of magnitude before saturating.

Figure 2(f) shows the characteristic curves $(t,S(t))$ corresponding to
different constant values of the tangent angle $\theta=\Theta^{(0)}$ [the
physical trajectories of these curves are shown in red in Fig. 2(a)]. Near the
tip ($S\approx 0$) characteristics initially converge towards one another,
implying a large range of tangent angles crowded into a small region, i.e.
sharpening. Previous discretizations of Eq. 2 interpreted this convergence as
a crossing of characteristics and thus the formation of a geometric shock. The
reformulated Eq. 4, however, reveals that characteristics bend away from one
another before ever crossing, thus preventing a finite-time blowup of
curvature. Characteristics farther from the tip are seen to change their
direction of travel, initially propagating outwards, and then inwards, before
they ultimately straighten and travel vertically. At late times, all
characteristics are seen to travel vertically, suggesting that a terminal
shape has arrived.

Exact solutions.— To examine the possibility of a terminal shape, we take a
$\theta$-derivative of Eq. 4 to obtain an evolution equation for the radius of
curvature [43, 10]:

$\frac{\partial R}{\partial
t}=V_{n}+\frac{\partial^{2}V_{n}}{\partial\theta^{2}}.$ (6)

Clearly, a steady-state of Eq. 6 is given by

$V_{n}=-V_{0}{}\sin{\theta}$ (7)

for any constant $V_{0}$, which is the recessional rate of the tip. Eq. 7
represents steady translation of a fixed shape. It is the only steady-state
$V_{n}$ that satisfies left-right symmetry. Inserting Eq. 7 into Eq. 5 and
inverting gives the equilibrium distribution of $R$,

$\frac{R^{*}}{R^{*}_{0}}=\frac{1+2\cos^{2}\theta}{\sin^{5}\theta}\,.$ (8)

This class of equilibrium solutions has one degree of freedom $R^{*}_{0}$,
which is the equilibrium radius of curvature at the tip. Exact expressions for
the Cartesian coordinates of this surface, along with solutions for the
corresponding axisymmetric (3D) problem, are given in the SI. Though
differences exist in the $\theta$–$L$ formulation of the 2D and 3D problems,
the final equilibrium solutions are identical when written in Cartesian
coordinates.

Figure 3:  Convergence towards the equilibrium morphology. (a) Left:
Overlaying interfaces in Fig. 2(a) shows a common shape to emerge at late
times. Right: Zooming-in near the apex further reveals the convergence towards
an equilibrium. (b) The rescaled radius of curvature $R/R_{0}$ tends to the
exact distribution predicted by Eq. 8. (c) Choosing initial shapes near the
equilibrium can lead to sharpening or blunting, as predicted by Eq. 10. Inset:
physical shape evolution of the case $\gamma=5/3$ reveals straight
characteristic paths.

To test the convergence to this final shape, Fig. 3(a) shows the simulated
interfaces from the previous example, but shifted to have the same apex. As
seen here and in the close-up, the interfaces indeed collapse to a single
profile at late times. Figure 3(b) shows that the corresponding distributions
of rescaled curvature-radius, $R(\theta,t)/R_{0}(t)$, converge to the
equilibrium shape Eq. 8.

Having observed the convergence to the predicted morphology, several questions
remain: What happens for different initial conditions? What determines the
final tip radius $R^{*}_{0}=\lim_{t\to\infty}R(0,t)$? And how can the results
be reconciled with previous infinite-curvature predictions [9, 26]? To address
these questions, we consider a local expansion in small $w=\cos\theta$:

$s(w,t)=a_{1}(t)w+a_{3}(t)w^{3}+\dots$ (9)

where odd-symmetry has been used. Thanks to the change of variables, $V_{n}$
can be calculated exactly for any power $w^{n}$. Inserting into Eq. 4
produces, at leading order, $\dot{R}_{0}\propto-R^{-1/4}_{0}$, which is
consistent with [9] and predicts finite-time blowup of curvature. However,
retaining the higher-order terms gives

$\dot{R}_{0}\propto-R^{-1/4}_{0}\left(1-\frac{3}{5}\gamma\right),\quad\gamma(t)=\frac{a_{3}(t)}{a_{1}(t)}\,,$
(10)

which is an exact relation (no truncation). Equation 10 opens the possibility
for the curvature divergence to be controlled by the term
$\left(1-\frac{3}{5}\gamma\right)$, and indeed the equilibrium solution Eq. 8
has the property $\gamma=5/3$.

To further examine this possibility, Fig. 3(c) shows the simulated dissolution
of three initial conditions surrounding the equilibrium:
$s(w,0)=a_{1}(0)w+a_{3}(0)w^{3}$, with $\gamma$ initially set to $5/3$, $4/3$,
and $3$. The figure confirms that $\gamma=5/3$ results in nearly constant
curvature 111The curvature does not remain exactly constant due to the fact
that this initial condition is not the terminal shape from Eq. 8, but a
truncation of it., whereas $\gamma>5/3$ (resp. $\gamma<5/3$) leads to
decreasing (resp. increasing) curvature, consistent with the sign of
$\dot{R}_{0}$ in Eq. 10. Thus, both tip sharpening and blunting are possible
[25, 26], with the value of $\gamma$ determining which occurs. The case
$\gamma=4/3$ leads to tip sharpening, but, due to the proximity to the
equilibrium, not nearly as much as in our first numerical example. Thus, the
enormous curvature growth observed in Fig. 2 should not always be expected, as
it depends on the initial shape.

We now turn attention to the experimentally-measured shapes that were shown in
Fig. 1, for planar (2D) and axisymmetric (3D) geometries. Figure 4(a) compares
the experimental profiles (shifted to have the same apex) to the equilibrium
morphology of Eq. 8 (thick gray curve). At late times, the experimental
profiles all collapse onto the predicted shape in both 2D and 3D, thus
conclusively confirming that Eq. 8 accurately describes the equilibrium spire-
morphology of a body dissolving under its own solute-induced convective flow.
This agreement with laboratory experiments also validates modeling assumptions
made, including the boundary-layer and quasi-steady approximations and the
omission of GT effects.

Figure 4:  Comparison with laboratory experiments shown in Fig. 1. (a) When
overlaid, the profiles measured from the planar (2D) and axisymmetric (3D)
experiments are seen to converge to the shape predicted by Eq. 8. Plotting the
experimentally measured surface coordinates on a log-scale confirms the far-
field prediction Eq. 11.

A second test is made possible by the far-field ($|s|\to\infty$) behavior of
the equilibrium solution

$\displaystyle y$
$\displaystyle\sim\frac{3}{4}\,{R_{0}^{*}}^{-1/3}\,x^{4/3}\,,$ (11)

which holds in both 2D and 3D [with $(x,y)\to(r,z)$ in 3D, see SI]. Figure
4(b) shows a log-scale comparison between this predicted 4/3-power law and the
experimental measurements. At late times, the experimental profiles indeed
converge to the predicted power law in both cases. We note that this power-law
is consistent with one of the similarity solutions found in [25, 26], which
would need to be asymptotically matched to an inner (near-tip) solution. Those
similarity solutions, however, predict continued evolution of shape, whereas
we have found convergence to a final form. Numerical and experimental evidence
suggests this final morphology to be a stable attractor.

Closer examination of our exact solutions offers an interpretation of the
flow-physics underlying the convergence in shape dynamics. Within the boundary
layer, a few competing effects exist. First, the apex is in contact with
nearly pure liquid, whereas a solute mixture washes over the downstream
portions. In isolation, this effect would cause the apex to retreat fastest.
On the other hand, the buoyancy-driven flow accelerates as it advances
downstream, due to the accumulation of dense solute as well as the increase in
surface steepness. This effect enhances convection-induced dissolution on
downstream portions. Which effect is stronger depends on the detailed geometry
of the object, and it is the interplay between the two that drives shape
change. Ultimately, balance is achieved by the steadily-translating
distribution, $V_{n}=-V_{0}{}\sin{\theta}$, which shows that the dissolution-
rate is highest at the tip ($\theta=\pi/2$) and decreases locally in
proportion to the surface steepness. At this stage, the mass loss rate of the
pinnacle has a simple scaling $dm/dt\sim
2\int_{\theta_{0}}^{\pi/2}V_{n}(\theta)R^{*}(\theta)\,d\theta\sim-
V_{0}R_{0}^{*}\sim-(R_{0}^{*})^{3/4}$, implying that the mass loss slows as
the tip sharpens.

Discussion.— In this Letter, we have described, in exact form, the final spire
morphology of a body being reshaped under its own dissolution-induced natural
convective flow, thereby concluding the search from [24, 38, 9, 25, 26].
Carefully designed numerics show that, rather than forming a geometric shock,
characteristics avoid crossing to pursue this terminal shape, which exhibits
large, but finite tip curvature. This situation is perhaps analogous to exact
solutions found in the context of free surface flows, whose finite curvature
reversed previous hypotheses on the formation of cusp singularities [45].

The simple, explicit nature of our solutions suggests that they may be used to
infer properties, e.g. age or past environmental conditions, of natural
structures. To take one example, suppose that a karst pinnacle at time $t_{0}$
has height $h(t_{0})$, width $d(t_{0})$, and that its apex dissolves at the
rate $\dot{h}(t_{0})=V(t_{0})$. As it nears the final shape, Eq. 11 suggests
$h(t)d(t)^{-4/3}=h(t_{0})d(t_{0})^{-4/3}$ and the constant tip velocity gives
$h(t)-h(t_{0})=V(t_{0})(t-t_{0})$. These two relationships comprise a closed
system for $(d(t),h(t))$ at any given time – including the past ($t<t_{0}$)
and the future ($t>t_{0}$)– thus offering the potential to estimate the past
dimensions or, if the dimensions can be estimated through other means, the age
of the structure. To take this idea one step further, the typical spacing $L$
between pinnacles in a stone forest approximates the initial width
$d(0)\approx L$, thus offering simple estimates for the pinnacle’s initial
height $h(0)=h(t_{0})[L/d(t_{0})]^{4/3}$ and its age
$t_{0}=[h(t_{0})-h(0)]/V(t_{0})$. Though natural systems involve a range of
other complicating factors (such as rainfall, turbulent boundary layers, and
fracture) our calculations, based principally on dissolution and fluid
dynamics, may offer a leading-order understanding of these amazing structures.

Aspects of our analysis can be extended to other physical systems. For
example, Eq. 6 can have a separable solution $R(\theta,t)=A(\theta)B(t)$,
corresponding to the self-similar evolution of erodible and soluble bodies
immersed in an externally forced flow [1, 42, 2]. Meanwhile, our approach can
be applied to dynamics with an opposite sign in $V_{n}$, seen in growing
systems like crystallization [43] and the formation of stalactites [17, 18].

## References

  * Ristroph  _et al._ [2012] L. Ristroph, M. N. J. Moore, S. Childress, M. J. Shelley, and J. Zhang, Sculpting of an erodible body by flowing water, Proc. Natl. Acad. Sci. 109, 19606 (2012).
  * Huang  _et al._ [2015] J. M. Huang, M. N. J. Moore, and L. Ristroph, Shape dynamics and scaling laws for a body dissolving in fluid flow, J. Fluid Mech. 765 (2015).
  * Sweeting [2012] M. M. Sweeting,  _Karst in China: its geomorphology and environment_ , Vol. 15 (Springer Science & Business Media, 2012).
  * Ford and Williams [2013] D. Ford and P. D. Williams,  _Karst hydrogeology and geomorphology_ (John Wiley & Sons, 2013).
  * Whitehead [1972] J. A. Whitehead, Moving heaters as a model of continental drift, Phys. Earth Planet. In. 5, 199 (1972).
  * Whitehead [1988] J. A. Whitehead, Fluid models of geological hotspots, Annu. Rev. Fluid Mech. 20, 61 (1988).
  * Zhong and Zhang [2005] J.-Q. Zhong and J. Zhang, Thermal convection with a freely moving top boundary, Phys. Fluids 17, 115105 (2005).
  * Huang  _et al._ [2018] J. M. Huang, J.-Q. Zhong, J. Zhang, and L. Mertz, Stochastic dynamics of fluid–structure interaction in turbulent thermal convection, J. Fluid Mech. 854 (2018).
  * Huang  _et al._ [2020] J. M. Huang, J. Tong, M. Shelley, and L. Ristroph, Ultra-sharp pinnacles sculpted by natural convective dissolution, Proc. Natl. Acad. Sci. 117, 23339 (2020).
  * Meakin and Jamtveit [2010] P. Meakin and B. Jamtveit, Geological pattern formation by growth and dissolution in aqueous systems, Proc. Math. Phys. Eng. Sci. 466, 659 (2010).
  * Dodds and Rothman [2000] P. S. Dodds and D. H. Rothman, Scaling, universality, and geomorphology, Annu. Rev. Earth Planet. Sci. 28, 571 (2000).
  * Short  _et al._ [2006] M. B. Short, J. C. Baygents, and R. E. Goldstein, A free-boundary theory for the shape of the ideal dripping icicle, Phys. Fluids 18, 083101 (2006).
  * Callander [1978] R. A. Callander, River meandering, Annu. Rev. Fluid Mech. 10, 129 (1978).
  * Stølum [1996] H.-H. Stølum, River meandering as self-organization process, Science 271, 1710 (1996).
  * Braudrick  _et al._ [2009] C. A. Braudrick, W. E. Dietrich, G. T. Leverich, and L. S. Sklar, Experimental evidence for the conditions necessary to sustain meandering in coarse-bedded rivers, Proc. Natl. Acad. Sci. USA 106, 16936 (2009).
  * Allison [1923] V. C. Allison, The growth of stalagmites and stalactites, J. Geol. 31, 106 (1923).
  * Short  _et al._ [2005a] M. B. Short, J. C. Baygents, J. W. Beck, D. A. Stone, R. S. Toomey III, and R. E. Goldstein, Stalactite growth as a free-boundary problem: A geometric law and its platonic ideal, Phys. Rev. Lett. 94, 018501 (2005a).
  * Short  _et al._ [2005b] M. B. Short, J. C. Baygents, and R. E. Goldstein, Stalactite growth as a free-boundary problem, Phys. Fluids 17, 083101 (2005b).
  * Amin  _et al._ [2019] K. Amin, J. Mac Huang, K. J. Hu, J. Zhang, and L. Ristroph, The role of shape-dependent flight stability in the origin of oriented meteorites, Proc. Natl. Acad. Sci. 116, 16180 (2019).
  * Ginés  _et al._ [2009] A. Ginés, M. Knez, T. Slabe, and W. Dreybrodt,  _Karst Rock Features. Karren sculpturing: Karren sculpturing_ , Vol. 9 (Založba ZRC, 2009).
  * Song [1986] L. H. Song, Origination of stone forests in China, Int. J. Speleol. 15, 1 (1986).
  * Veress  _et al._ [2008] M. Veress, D. Lóczy, Z. Zentai, G. Tóth, and R. Schläffer, The origin of the Bemaraha tsingy (Madagascar), Int. J. Speleol. 37, 6 (2008).
  * Slabe  _et al._ [2016] T. Slabe, A. Hada, and M. Knez, Laboratory modeling of karst phenomena and their rock relief on plaster: Subsoil karren, rain flutes karren and caves, Acta Carsologica 45 (2016).
  * Nakouzi  _et al._ [2014] E. Nakouzi, R. E. Goldstein, and O. Steinbock, Do dissolving objects converge to a universal shape?, Langmuir (2014).
  * Pegler and Wykes [2020] S. S. Pegler and M. S. D. Wykes, Shaping of melting and dissolving solids under natural convection, J. Fluid Mech. 900 (2020).
  * Pegler and Wykes [2021] S. S. Pegler and M. S. D. Wykes, The convective Stefan problem: shaping under natural convection, J. Fluid Mech. 915 (2021).
  * Aaron and Kotler [1970] H. B. Aaron and G. R. Kotler, The effects of curvature on the dissolution kinetics of spherical precipitates, Met. Sci. J. 4, 222 (1970).
  * Perez [2005] M. Perez, Gibbs–Thomson effects in phase transformations, Scr. Mater. 52, 709 (2005).
  * Moore [2017] M. N. J. Moore, Riemann-Hilbert Problems for the Shapes Formed by Bodies Dissolving, Melting, and Eroding in Fluid Flows, Comm. Pure Appl. Math. 70, 1810 (2017).
  * Huang  _et al._ [2021] J. M. Huang, M. J. Shelley, and D. B. Stein, A stable and accurate scheme for solving the Stefan problem coupled with natural convection using the Immersed Boundary Smooth Extension method, J. Comput. Phys. 432, 110162 (2021).
  * Rycroft and Bazant [2016] C. H. Rycroft and M. Z. Bazant, Asymmetric collapse by dissolution or melting in a uniform flow, Proc. Math. Phys. Eng. Sci. 472, 20150531 (2016).
  * Mitchell and Spagnolie [2017] W. H. Mitchell and S. E. Spagnolie, A generalized traction integral equation for Stokes flow, with applications to near-wall particle mobility and viscous erosion, J. Comput. Phys. 333, 462 (2017).
  * Hewett and Sellier [2017] J. N. Hewett and M. Sellier, Evolution of an eroding cylinder in single and lattice arrangements, J. Fluids Struct. 70, 295 (2017).
  * Derr  _et al._ [2020] N. J. Derr, D. C. Fronk, C. A. Weber, A. Mahadevan, C. H. Rycroft, and L. Mahadevan, Flow-driven branching in a frangible porous medium, Phys. Rev. Lett. 125, 158002 (2020).
  * Quaife and Moore [2018] B. D. Quaife and M. N. J. Moore, A boundary-integral framework to simulate viscous erosion of a porous medium, J. Comput. Phys. 375, 1 (2018).
  * Chiu  _et al._ [2020] S.-H. Chiu, M. N. J. Moore, and B. Quaife, Viscous transport in eroding porous media, J. Fluid Mech. 893 (2020).
  * Ladd  _et al._ [2020] A. J. C. Ladd, L. Yu, and P. Szymczak, Dissolution of a cylindrical disk in Hele-Shaw flow: a conformal-mapping approach, J. Fluid Mech. 903 (2020).
  * Wykes  _et al._ [2018] M. S. D. Wykes, J. M. Huang, G. A. Hajjar, and L. Ristroph, Self-sculpting of a dissolvable body due to gravitational convection, Phys. Rev. Fluids 3, 043801 (2018).
  * Schlichting and Gersten [2016] H. Schlichting and K. Gersten,  _Boundary-Layer Theory_ (Springer, 2016).
  * Hou  _et al._ [1994] T. Y. Hou, J. S. Lowengrub, and M. J. Shelley, Removing the stiffness from interfacial flows with surface tension, J. Comput. Phys. 114, 312 (1994).
  * Alben  _et al._ [2002] S. Alben, M. Shelley, and J. Zhang, Drag reduction through self-similar bending of a flexible body, Nature 420, 479 (2002).
  * Moore  _et al._ [2013] M. N. J. Moore, L. Ristroph, S. Childress, J. Zhang, and M. J. Shelley, Self-similar evolution of a body eroding in a fluid flow, Phys. Fluids 25, 116602 (2013).
  * Wettlaufer  _et al._ [1994] J. S. Wettlaufer, M. Jackson, and M. Elbaum, A geometric model for anisotropic crystal growth, J. Phys. A: Math. Gen. 27, 5957 (1994).
  * Note [1] The curvature does not remain exactly constant due to the fact that this initial condition is not the terminal shape from Eq. 8, but a truncation of it.
  * Jeong and Moffatt [1992] J.-T. Jeong and H. Moffatt, Free-surface cusps associated with flow at low Reynolds number, J. Fluid Mech. 241, 1 (1992).