# Residual-Concatenate Neural Network with Deep Regularization Layers for
Binary Classification

Abhishek Gupta
Research Scholar
University of Mumbai
Mumbai, MH, India
<EMAIL_ADDRESS>
Sruthi Nair
Master of Engineering
Vidyalankar Institute of Technology
Mumbai, MH, India
<EMAIL_ADDRESS>
Raunak Joshi
Mentor
University of Mumbai
Mumbai, MH, India
<EMAIL_ADDRESS>
Vidya Chitre
Assistant Professor - Department of IT
Vidyalankar Institute of Technology
Mumbai, MH, India
<EMAIL_ADDRESS>

###### Abstract

Many complex Deep Learning models are used with different variations for
various prognostication tasks. The higher learning parameters not necessarily
ensure great accuracy. This can be solved by considering changes in very deep
models with many regularization based techniques. In this paper we train a
deep neural network that uses many regularization layers with residual and
concatenation process for best fit with Polycystic Ovary Syndrome Diagnosis
prognostication. The network was built with improvements from every step of
failure to meet the needs of the data and achieves an accuracy of 99.3%
seamlessly.

_Keywords_ Concatenation $\cdot$ Regularization $\cdot$ Residuals

##  1 Introduction

Deep Learning(LeCun et al., 2015) definitely has made many advancements in
last decade by introducing many different types of neural networks. These
include solving different problems from the domain of vision, text and audio.
The area of deep learning has been using neural networks which is a learning
procedure that replicates the biological human neuron system. It learns
representations and patters from the data like humans do. These
representations are highly adaptable with changes in the learning patterns to
best fit the needs of the problem. The standard neural network has an input
layer and an output layer with $n$ number of hidden layers. These hidden
layers are responsible for the learning process that takes place and gives to
the output layer. The entire process of any neural network is a two-step
process consisting of forward propagation(Hirasawa et al., 1996) and backward
propagation(Rumelhart et al., 1986). The forward propagation is responsible
for learning patterns from the input data and backward propagation is
responsible for improving the losses incurred in the learning process of
forward propagation. The input layer considers the data in form of tensor and
initializes the weight and bias for the network. These parameters are passed
later to the hidden layers. These hidden layers are responsible for modulating
these weight and biases with respect to input tensor and helps learn
representations. These weight and biases alone are not effective without the
use of activation functions(Nwankpa et al., 2018) which are a set of
mathematical functions that purely focus on the aspect of considering the
right aspect of information from the entire representation. These activation
functions differ from task to task, considering binary classification(Kumari
and Srivastava, 2017), an activation function like sigmoid(Minai and Williams,
1993) will be used that differentiates the 2 classes very effectively.
Assuming classes in the labels are more than 2, softmax(Bridle, 1989)(Uffink,
2017) is used which gives probabilities of all the predictable classes where
the class with the highest probability can be extracted very efficiently. For
learning representations from the set of hidden layers that are interconnected
an activation function like tanh(Ryck et al., 2021) or non-linearity
activation like rectified linear unit also known as ReLU(Agarap, 2018) is
used. Amidst this learning process the loss is calculated with a loss based
function. These loss functions compute and keep track of the information lost
during learning representations. This completes the process of forward
propagation and process is transferred to backward propagation where the
gradients for loss are computed and changes in learnt weights and biases are
computed. This entire process consists of one epoch(Afaq and Rao, 2020). Now
in this paper we focus particularly on binary classification task with a deep
neural network. When the number of hidden layers are very high in number the
neural network is considered as deep neural network. Now training a very deep
neural network can run into a problem of learning representations either more
or less. This is where regularization techniques are introduced and these
prove to be a better for improving the performance. Considering the binary
classification problem, we had a wide horizon for selecting data. In this case
we considered working with Polycystic Ovary Syndrome diagnosis dataset. This
is disorder observed in women and the data is very high dimensional and
conducive for binary classification. Machine Learning is widely used for PCOS
prognostication in many different contexts. Considering most commonly used
algorithms of linear models or distance based algorithms(Chauhan et al., 2021)
has been done. Apart from that work with Discriminant Analysis(Gupta et al.,
2022) is done. Extending to ensemble learning(Denny et al., 2019) and the
varieties it provides, the work with bagging(Kanvinde et al., 2022),
boosting(Gupta et al., 2021) and stacking(Nair et al., 2022) has been
performed effectively. In this paper we specifically target the weaker
sections of all state-of-the-art approaches and try to come up with techniques
to get past the transgressions of learning methods.

##  2 Methodology

###  2.1 Basic Network

The basic network is build to have a base to start with for making changes.
The data used has 41 dimensions of feature set and 1-dimensional vector of
labels. Training of 3-5 hidden layers is done specifically that targets the
basic learning with one output layer with 1 hidden neuron and sigmoid
activation function. The hidden neurons from hidden layers are highly
arbitrary and loss function used is binary cross-entropy with Adam(Kingma and
Ba, 2017) loss optimization(Reddy and Reddy, ) function. Now this network is
particularly not efficient for learning representations and runs into an
under-fitting(Zhang et al., 2019) problem. So this is where a more deeper
network needs to be trained so the number of hidden layers will be increased
but the model will give high variance(Mehta et al., 2019) and run into an
over-fitting(Zhang et al., 2019) problem. For this the need of regularization
arises.

###  2.2 Weight Decay

The weight decay(Xie et al., 2021) process is also known as regularization and
there are 2 types, $L_{1}$ and $L_{2}$(Ng, 2004)(Szabo and Lorincz, 2004)
regularization and it directly works for solving the over-fitting problem of
the networks. Considering the $L_{2}$ regularization first, the regularization
parameter is Euclidean Norm(Wu et al., 2019) and can be defined by the sign
$\Omega$. It considers the weight matrix of the learning representation and
computes sum over all the squared weight values of it. The trained network can
be seen in the Fig. 1

Figure 1: L2 Regularization with Batch Normalization

The batch normalization layer has non-learnable parameters that cuts off the
excess slack for any neural network learning process very efficiently. The
formula for the $L_{2}$ Regularization is given as

$\Omega(W)=\sum_{i}\sum_{j}W^{2}_{ij}$ (1)

where the W is the weight and now this further changes the loss function. One
scalar parameter of $\alpha$ is used. This changes the loss function and can
be given as

$L_{2}(W)=\frac{\alpha}{2}\sum_{i}\sum_{j}W^{2}_{ij}+L(W)$ (2)

where the $L_{2}$ is regularized loss function and along with the Euclidean
Norm the loss function is added and changes are made. Similarly same
architecture with $L_{1}$ Regularization is used. It considers the absolute
values of the weight parameters. The formula for this can be given as

$\Omega(W)=\sum_{i}\sum_{j}|W_{ij}|$

and loss function can be updated as

$L_{1}(W)=\alpha\sum_{i}\sum_{j}|W_{ij}|+L(W)$ (3)

Now these architectures are not performing efficiently and the changes can be
made using a different procedure.

###  2.3 Concatenation

This is a very different process where you combine best of both networks. In
this the network are trained in a parallel style. Assume training multiple
networks simultaneously and the outputs are received at other end where they
are combined and given to the output layer. This has an advantage as well as
disadvantage. The advantage is that results can be improved as more parameters
are learnt by network. Disadvantage can be that if one network is not able to
perform effectively, the results can be impacted. But this procedure proves to
be effective in most of the cases.

Figure 2: Regularization Model with Concatenation

The Fig. 2 is the graphical representation of the model. The input
simultaneously imparts the data to 2 networks independently which train in the
parallel fashion and give their respective outputs to the concatenation(Noreen
et al., 2020) layer. The concatenation layer has the rule that the dimensions
of both the networks are supposed to be same. We can consider the number of
hidden neurons in layman terms. If they are not same the operation is not
performed effectively. Later the output from the concatenation layer is given
to the output layer. This process increases the number of trainable
parameters. With this type of network we are able to curb the under-fitting
problem but now the network runs into an over-fitting problem. Since the
network trained is very deep and fits on the data very efficiently and does
not consider the validation data effectively, more new variation needs to be
introduced.

###  2.4 Residuals

Right before understanding the concept of Residual, the final architecture can
be represented with Fig. 3

Figure 3: Residual Concatenated Model

The residual(He et al., 2016) blocks are the main essence of this
architecture. The over-fitting issue persists because of the saturation in the
learning process. Deeper the layers, more representations will be computed and
the calculating the gradients in the back propagation becomes more complex.
The residual block inspired by ResNet(He et al., 2016), works with the concept
of skip connections or residual connections. The skip connections(Wu et al.,
2020) particularly work when the degradation problem starts. The batch
normalization(Ioffe and Szegedy, 2015) layer was eliminated for the fact that
network was over-fitting. Highway networks(Srivastava et al., 2015) was the
concept from which the residual blocks were derived. The calculable parameters
are important where the $L_{1}$ model has 809472 parameters, whereas $L_{2}$
model has the same parameters. The total amount of parameters of combined
network are 1750273 parameters. The activation functions used in the hidden
layer are ReLU(Agarap, 2018). The number of hidden neurons are 512 in every
single layer. The output from both the regularization models is applied with a
residual block and the output is applied with a ReLU layer. Later the output
from ReLU layer is given to the concatenation layer. The output of the
concatenation layer is given to a 128 hidden neuron layer with ReLU
activation. This finally gives output to the last layer with sigmoid
activation function.

##  3 Results

The main inferences can be derived from the result section. The results
consider results of all the different models used and what was the need for
using a better model.

Figure 4: Training Accuracies

The accuracy and loss can be visualized to have an intuition of how the
performance of the model is taking place. Similarly the graph for loss can
also be visualized that shows the loss reduced with respect to number of
epochs with all the model variations.

Figure 5: Training Losses

The Fig. 4 and 5 give very distinct visual representations of the model
training performance, though the question arises how come we assume the
residual model performs better even though the training accuracy is lower as
compared to other models. The answer to this is the differences between
training and validation accuracy. It is a key determining factor for defining
the state of the model. If the training accuracy is way lower than expected,
it is said to be in state of under-fitting, while if accuracy of training is
way higher and testing accuracy is pretty much stagnant, inference can be
drawn that model is in state of over-fitting. The accuracies of the other
models is 100% and convergence of the loss is not as good as expected. This is
the need for the residual model.

Table 1: Training and Validation Accuracies Algorithm | Training Accuracy | Validation Accuracy
---|---|---
$L_{1}$-Reg | 100% | 92.72%
$L_{2}$-Reg | 100% | 94.54%
Concat | 100% | 92.72%
Residual | 99.3% | 95.4%

The Table 1 gives the proof for the need of the Residual changes. The other
models are over-fitting and residual block network did performed better than
the others. This pretty much eliminated the need of observing the values of
the loss table. But just for the sake of deriving the perfect inference, let
us look at Table 2.

Table 2: Training and Validation Losses Algorithm | Training Loss | Validation Loss
---|---|---
$L_{1}$-Reg | 11.91% | 28.81%
$L_{2}$-Reg | 11.73% | 25.9%
Concat | 17.65% | 36.13%
Residual | 8.11% | 21.48%

When we consider the losses, residual performs the best amongst others. The
need for residual was mandatory based on the performance metrics of other
variations of neural network.

##  4 Conclusion

The major advancements in the area of deep learning were referred to create a
very subtle architecture that points out some transgressions with binary
classification problem. The demonstration of the approach was done starting
from very basic of the neural networks. The data has 41 dimensional feature
set. In the start a very simple 3-5 layered architecture was designed which
was increased and regularization was included. $L_{1}$ and $L_{2}$
regularization was introduced that has some issues which were later addressed
with concatenation model. This model has high number of parameters as it is
combined to give the best of 2 worlds. Yet some over-fitting issues were
experienced for which residual based blocks were introduced that gave the best
performance. This paper targets a very subtle implementation of neural
networks for high dimensional binary classification problem. In future by
referring this paper, many new methods can be formulated for which we would be
more than happy.

## Acknowledgment

Our team thanks Mr. Prasoon Kottarathil for his effort in gathering and
providing the dataset of Polycystic Ovary Syndrome on Kaggle.

## References

  * LeCun et al. [2015] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  _nature_ , 521(7553):436–444, 2015.
  * Hirasawa et al. [1996] K. Hirasawa, M. Ohbayashi, M. Koga, and M. Harada.  Forward propagation universal learning network.  In _Proceedings of International Conference on Neural Networks (ICNN’96)_ , volume 1, pages 353–358 vol.1, 1996.  doi:10.1109/ICNN.1996.548917.
  * Rumelhart et al. [1986] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning representations by back-propagating errors.  _Nature_ , 323:533–536, 1986.
  * Nwankpa et al. [2018] Chigozie Nwankpa, Winifred Ijomah, Anthony Gachagan, and Stephen Marshall.  Activation functions: Comparison of trends in practice and research for deep learning, 2018.
  * Kumari and Srivastava [2017] Roshan Kumari and Saurabh Kr. Srivastava.  Machine learning: A review on binary classification.  _International Journal of Computer Applications_ , 160(7):11–15, Feb 2017.  ISSN 0975-8887.  doi:10.5120/ijca2017913083.  URL http://www.ijcaonline.org/archives/volume160/number7/27084-2017913083.
  * Minai and Williams [1993] Ali A. Minai and Ronald D. Williams.  Original contribution: On the derivatives of the sigmoid.  _Neural Netw._ , 6(6):845–853, jun 1993.  ISSN 0893-6080.  doi:10.1016/S0893-6080(05)80129-7.  URL https://doi.org/10.1016/S0893-6080(05)80129-7.
  * Bridle [1989] John Bridle.  Training stochastic model recognition algorithms as networks can lead to maximum mutual information estimation of parameters.  In D. Touretzky, editor, _Advances in Neural Information Processing Systems_ , volume 2. Morgan-Kaufmann, 1989.
  * Uffink [2017] Jos Uffink.  Boltzmann’s Work in Statistical Physics.  In Edward N. Zalta, editor, _The Stanford Encyclopedia of Philosophy_. Metaphysics Research Lab, Stanford University, Spring 2017 edition, 2017.
  * Ryck et al. [2021] Tim De Ryck, Samuel Lanthaler, and Siddhartha Mishra.  On the approximation of functions by tanh neural networks.  _Neural networks : the official journal of the International Neural Network Society_ , 143:732–750, 2021.
  * Agarap [2018] Abien Fred Agarap.  Deep learning using rectified linear units (relu).  _arXiv preprint arXiv:1803.08375_ , 2018.
  * Afaq and Rao [2020] Saahil Afaq and Smitha Rao.  Significance of epochs on training a neural network.  _International Journal of Scientific & Technology Research_, 9:485–488, 2020.
  * Chauhan et al. [2021] Preeti Chauhan, Pooja Patil, Neha Rane, Pooja Raundale, and Harshil Kanakia.  Comparative analysis of machine learning algorithms for prediction of pcos.  In _2021 International Conference on Communication information and Computing Technology (ICCICT)_ , pages 1–7, 2021.  doi:10.1109/ICCICT50803.2021.9510128.
  * Gupta et al. [2022] Abhishek Gupta, Himanshu Soni, Raunak Joshi, and Ronald Melwin Laban.  Discriminant analysis in contrasting dimensions for polycystic ovary syndrome prognostication.  _arXiv preprint arXiv:2201.03029_ , 2022.
  * Denny et al. [2019] Amsy Denny, Anita Raj, Ashi Ashok, C Maneesh Ram, and Remya George.  i-hope: Detection and prediction system for polycystic ovary syndrome (pcos) using machine learning techniques.  In _TENCON 2019 - 2019 IEEE Region 10 Conference (TENCON)_ , pages 673–678, 2019.  doi:10.1109/TENCON.2019.8929674.
  * Kanvinde et al. [2022] Nandan Kanvinde, Abhishek Gupta, and Raunak Joshi.  Binary classification for high dimensional data using supervised non-parametric ensemble method.  _arXiv preprint arXiv:2202.07779_ , 2022.
  * Gupta et al. [2021] Abhishek M Gupta, Sannidhi S Shetty, Raunak M Joshi, and Ronald Melwin Laban.  Succinct differentiation of disparate boosting ensemble learning methods for prognostication of polycystic ovary syndrome diagnosis.  In _2021 International Conference on Advances in Computing, Communication, and Control (ICAC3)_ , pages 1–5. IEEE, 2021.
  * Nair et al. [2022] Sruthi Nair, Abhishek Gupta, Raunak Joshi, and Vidya Chitre.  Combining varied learners for binary classification using stacked generalization.  _arXiv preprint arXiv:2202.08910_ , 2022.
  * Kingma and Ba [2017] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * [19] SVG Reddy and K Thammi Reddy.  Optimization of deep learning using various optimizers, loss functions and dropout.
  * Zhang et al. [2019] Haotian Zhang, Lin Zhang, and Yuan Jiang.  Overfitting and underfitting analysis for deep learning based end-to-end communication systems.  In _2019 11th International Conference on Wireless Communications and Signal Processing (WCSP)_ , pages 1–6, 2019.  doi:10.1109/WCSP.2019.8927876.
  * Mehta et al. [2019] Pankaj Mehta, Marin Bukov, Ching-Hao Wang, Alexandre G. R. Day, Clint Richardson, Charles K. Fisher, and David J. Schwab.  A high-bias, low-variance introduction to machine learning for physicists.  _Physics reports_ , 810:1–124, 2019.
  * Xie et al. [2021] Zeke Xie, Issei Sato, and Masashi Sugiyama.  Understanding and scheduling weight decay, 2021.
  * Ng [2004] A. Ng.  Feature selection, l1 vs. l2 regularization, and rotational invariance.  _Proceedings of the twenty-first international conference on Machine learning_ , 2004.
  * Szabo and Lorincz [2004] Z. Szabo and A. Lorincz.  L1 regularization is better than l2 for learning and predicting chaotic systems, 2004.
  * Wu et al. [2019] Yan Wu, Mihaela Rosca, and Timothy Lillicrap.  Deep compressed sensing.  In _International Conference on Machine Learning_ , pages 6850–6860. PMLR, 2019.
  * Noreen et al. [2020] Neelum Noreen, Sellappan Palaniappan, Abdul Qayyum, Iftikhar Ahmad, Muhammad Imran, and Muhammad Shoaib.  A deep learning model based on concatenation approach for the diagnosis of brain tumor.  _IEEE Access_ , 8:55135–55144, 2020.
  * He et al. [2016] 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.
  * Wu et al. [2020] Dongxian Wu, Yisen Wang, Shu-Tao Xia, James Bailey, and Xingjun Ma.  Skip connections matter: On the transferability of adversarial examples generated with resnets, 2020.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International conference on machine learning_ , pages 448–456. PMLR, 2015.
  * Srivastava et al. [2015] Rupesh Kumar Srivastava, Klaus Greff, and Jürgen Schmidhuber.  Highway networks.  _arXiv preprint arXiv:1505.00387_ , 2015.


# Reply to “Comment on: ‘Case for a U(1)π Quantum Spin Liquid Ground State
in the Dipole-Octupole Pyrochlore Ce2Zr2O7’ ”

E. M. Smith  Department of Physics and Astronomy, McMaster University,
Hamilton, Ontario L8S 4M1, Canada  Brockhouse Institute for Materials
Research, McMaster University, Hamilton, Ontario L8S 4M1, Canada  O. Benton
Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Straße 38,
Dresden 01187, Germany  D. R. Yahne  Department of Physics, Colorado State
University, 200 West Lake Street, Fort Collins, Colorado 80523-1875, USA  B.
Placke  Max Planck Institute for the Physics of Complex Systems, Nöthnitzer
Straße 38, Dresden 01187, Germany  R. Schäfer  Max Planck Institute for the
Physics of Complex Systems, Nöthnitzer Straße 38, Dresden 01187, Germany  J.
Gaudet  Department of Physics and Astronomy, McMaster University, Hamilton,
Ontario L8S 4M1, Canada  J. Dudemaine  Département de Physique, Université de
Montréal, Montréal, Canada  Regroupement Québécois sur les Matériaux de Pointe
(RQMP)  A. Fitterman  Département de Physique, Université de Montréal,
Montréal, Canada  Regroupement Québécois sur les Matériaux de Pointe (RQMP)
J. Beare  Department of Physics and Astronomy, McMaster University, Hamilton,
Ontario L8S 4M1, Canada  A. R. Wildes  Institut Laue-Langevin, 71 Avenue des
Martyrs CS 20156, 38042 Grenoble Cedex 9, France  S. Bhattacharya  Université
Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405 Orsay, France
T. DeLazzer  Department of Physics, Colorado State University, 200 West Lake
Street, Fort Collins, Colorado 80523-1875, USA  C. R. C. Buhariwalla
Department of Physics and Astronomy, McMaster University, Hamilton, Ontario
L8S 4M1, Canada  N. P. Butch  Center for Neutron Research, National Institute
of Standards and Technology, MS 6100 Gaithersburg, Maryland 20899, USA  R.
Movshovich  Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
J. D. Garrett  Brockhouse Institute for Materials Research, McMaster
University, Hamilton, Ontario L8S 4M1, Canada  C. A. Marjerrison  Brockhouse
Institute for Materials Research, McMaster University, Hamilton, Ontario L8S
4M1, Canada  J. P. Clancy  Department of Physics and Astronomy, McMaster
University, Hamilton, Ontario L8S 4M1, Canada  Brockhouse Institute for
Materials Research, McMaster University, Hamilton, Ontario L8S 4M1, Canada  E.
Kermarrec  Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides,
91405 Orsay, France  G. M. Luke  Department of Physics and Astronomy, McMaster
University, Hamilton, Ontario L8S 4M1, Canada  Brockhouse Institute for
Materials Research, McMaster University, Hamilton, Ontario L8S 4M1, Canada  A.
D. Bianchi  Département de Physique, Université de Montréal, Montréal, Canada
Regroupement Québécois sur les Matériaux de Pointe (RQMP)  K. A. Ross
Department of Physics, Colorado State University, 200 West Lake Street, Fort
Collins, Colorado 80523-1875, USA  Canadian Institute for Advanced Research,
661 University Avenue, Toronto, Ontario M5G 1M1, Canada.  B. D. Gaulin
Department of Physics and Astronomy, McMaster University, Hamilton, Ontario
L8S 4M1, Canada  Brockhouse Institute for Materials Research, McMaster
University, Hamilton, Ontario L8S 4M1, Canada  Canadian Institute for Advanced
Research, 661 University Avenue, Toronto, Ontario M5G 1M1, Canada.

###### Abstract

In his comment [arXiv:2209:03235], S. W. Lovesey argues that our analysis of
neutron scattering experiments performed on Ce2Zr2O7 is invalid. Lovesey
argues that we have not properly accounted for the higher-order multipolar
contributions to the magnetic scattering and that our use of pseudospin-$1/2$
operators to describe the scattering is inappropriate. In this reply, we show
that the multipolar corrections discussed by Lovesey only become significant
at scattering wavevectors exceeding those accessed in our experiments. This in
no way contradicts or undermines our work, which never claimed a direct
observation of scattering from higher-order multipoles. We further show that
Lovesey’s objections to our use of pseudospins are unfounded, and that the
pseudospin operators are able to describe all magnetic scattering processes at
the energy scale of our experiments, far below the crystal field gap. Finally,
we comment on certain assumptions in Lovesey’s calculations of the scattering
amplitude which are inconsistent with experiment.

The comment [1] presents two principal objections to our original study [2]:

  1. 1.

It is argued that we did not properly account for corrections to the neutron
scattering cross section coming from higher-order multipoles (beyond the
dipolar approximation).

  2. 2.

It is argued that using pseudospin degrees of freedom to describe magnetic
scattering is fundamentally inappropriate, since the pseudospins “are
irrelevant in the material world” [1].

In this response, we will show that:

  1. 1.

Higher order multipolar corrections to the scattering intensity are
insignificant over the range of scattering wavevectors probed in our
experiments. This can be seen by calculating the form factors for the higher
order contributions [see Fig. 1]. Directly observing the higher multipolar
scattering was not an aim of our study, and we did not claim any such
observation. As such, our analysis of the scattering, in which we work within
the dipole approximation, is appropriate.

  2. 2.

The use of pseudospin-$1/2$ operators, which is standard in the literature on
rare-earth frustrated magnetism [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] , is the
appropriate way to treat the low energy magnetic degrees of freedom in
Ce2Zr2O7. Far from being “irrelevant”, these operators are directly related to
measurable quantities and can account for all magnetic scattering processes
(including multipolar scattering) which involve only the states of the low
energy Kramers doublet (which is all of the magnetic scattering presented in
[2]).

  3. 3.

The calculation of the scattering amplitude in [1] effectively assumes time-
reversal symmetry breaking order, which is inconsistent with the experimental
results.

##  I Goals and outcomes of our study

It appears to us that the comment [1] is partly motivated by a
misunderstanding of the claims in [2]. We therefore begin the response by
briefly reviewing what we actually set out to do, and achieved, in [2]. This
will then set the context for the remainder of our response.

Ce2Zr2O7 was identified in previous experimental studies as a highly
frustrated magnet and quantum spin liquid (QSL) candidate [14, 15]. Neutron
scattering showed that the crystal electric fields (CEF) on the Ce3+ sites
separate a low energy Kramers doublet from the rest of the CEF spectrum, by a
large gap $\Delta_{\rm CEF}\approx 56{\rm meV}$ [14, 15]. Characterization of
the wave-functions of the low energy CEF states revealed that they transform
according to the $\Gamma_{5}^{+}$ and $\Gamma_{6}^{+}$ irreducible
representations of the $D_{3d}$ point group. This characterisation identifies
Ce2Zr2O7 as belonging to the family of dipolar-octupolar pyrochlores, whose
novel properties were originally predicted in [16].

In particular, the symmetry allowed nearest-neighbor exchange interactions in
such a material lead to a Hamiltonian [Eq. (1) of [2]] which has been
predicted to give rise to a series of QSL states, occupying a large proportion
of parameter space [17, 18]. Given this context, the aims of our study were:

  1. 1.

To further characterize the low temperature state of Ce2Zr2O7 using heat
capacity, susceptibility and neutron scattering measurements

  2. 2.

To use these measurements to constrain the exchange parameters entering into
the Hamiltonian describing Ce2Zr2O7.

  3. 3.

To place those estimates within published phase diagrams of the Hamiltonian,
and thereby arrive at an inference of the ground state.

These aims were achieved, resulting in the prediction that the ground state of
Ce2Zr2O7 is a $U(1)_{\pi}$ spin liquid.

We emphasise that, contrary to the depiction in [1], it was not our goal to
provide further evidence that Ce2Zr2O7 is a dipolar-octupolar pyrochlore. We
regard this as having already been proven by the experiments measuring the CEF
states [14, 15], and we simply assume it for the purpose of our study.

Further, again contrary to [1], it was also not our goal to directly measure
the higher order multipolar contribution to the neutron scattering intensity
(unlike e.g. [11]) and we did not claim to do so. In order to do this, we
would have needed to measure scattering at much larger scattering wavevectors
(see Section II). Measuring such scattering was not necessary for the purposes
for which we used the neutron scattering data in our study, as discussed
below.

We now turn to discuss the role that the theoretical calculations of the
neutron scattering response played in our study. The initial estimation of the
exchange parameters was actually made by fitting heat capacity and
susceptibility measurements to model calculations, not by using the scattering
data (again, contrary to the description of our work in [1]). It was only once
optimal sets of parameters had already been determined, that we compared model
calculations of the neutron scattering intensity to measurements.

This comparison served two purposes: to differentiate between two local optima
in the goodness-of-fit with respect to the thermodynamic measurements and as
an additional check on the consistency of our theory. For the first purpose,
it was sufficient to use data at relatively small momentum transfer (in the
range $[0.43,0.96]{\rm\AA}^{-1}$, see Fig. 10 of [2]). For the second purpose,
we used all the data at hand, and this involved a range of momenta $q\lesssim
2.5{\rm\AA}^{-1}$. Although at such wavevectors the scattering is only
meaningfully sensitive to magnetic dipoles (see Section II) this is
nevertheless a nontrivial test of the theory as the wavevector and energy
dependence of the dipole correlations depends on all exchange parameters. It
would be interesting to test the theory with data from larger scattering
wavevectors, but that is beyond the scope of [2].

Figure 1: Squares of form factors multiplying different multipolar
contributions to the scattering amplitude for magnetic scattering from a Ce3+
ion, over the momentum range $q\in[0,2.5]{\rm\AA}^{-1}$ relevant to the
experiments in [2]. The form factor arising from the dipolar approximation
$f_{\rm dip}(q)$ [Eq. (2)], used in our calculation, agrees closely with
$f^{\prime}(q)$, the form factor of the leading, dipolar, contribution in Eq.
(5). The form factors of the higher order multipoles $h(q)$ and $g(q)$ are
both small over this momentum range. This justifies the use of the dipolar
approximation in theoretical calculations of the scattering, and demonstrates
that the higher order multipole corrections discussed in [1] will become
significant only at larger wavevectors. Calculations are based on analytical
approximations to the functions $\langle j_{n}(q)\rangle$ given in [19].

##  II Negligibility of higher-order multipolar corrections to the scattering
cross-section

We now consider the neutron scattering cross section and show that the
multipolar corrections discussed in [1] are negligibly small over the range of
wavevectors we access.

The magnetic scattering cross section is related to the correlations of the
Fourier transform of the magnetization operator $\hat{{\bf M}}({\bf q})$:

$\displaystyle\frac{d\sigma}{d\Omega}\propto\sum_{\alpha\beta}\left(\delta_{\alpha\beta}-\frac{q_{\alpha}q_{\beta}}{q^{2}}\right)\langle\hat{M}_{\alpha}(-{\bf
q})\hat{M}_{\beta}({\bf q})\rangle$ (1)

where $\left(\delta_{\alpha\beta}-\frac{q_{\alpha}q_{\beta}}{q^{2}}\right)$
projects out the longitudinal component of $\hat{{\bf M}}$. When the
wavevector ${\bf q}$ is small enough that the scattering is not sensitive to
atomic scale variations of the magnetisation, we can rewrite $\hat{{\bf
M}}({\bf q})$ as:

$\displaystyle\hat{{\bf M}}({\bf q})\approx f_{\rm dip}(q){\bf m}({\bf q})$
(2)

where ${\bf m}({\bf q})$ is the lattice Fourier transform of the ionic
magnetic moments and $f_{\rm dip}(q)$ is a scalar form factor. This is the
dipole approximation [20, 21]. In this approximation the form factor is:

$f_{\rm dip}(q)=\langle j_{0}(q)\rangle+\frac{(2-g_{J})}{g_{J}}\langle
j_{2}(q)\rangle$ (3)

where $g_{J}$ is the Landé $g$-factor ($=\frac{6}{7}$ for Ce3+) and $\langle
j_{n}(q)\rangle$ are the expectation values of spherical Bessel functions
$J_{n}(q\ r)$ evaluated with respect to the radial electron density of the
magnetic ions.

This leads to the following result for the scattering cross-section, from
which Eq. (F1) of [2] can be derived (see Section III):

$\displaystyle\frac{d\sigma}{d\Omega}\propto|f_{\rm
dip}(q)|^{2}\sum_{\alpha\beta}\left(\delta_{\alpha\beta}-\frac{q_{\alpha}q_{\beta}}{q^{2}}\right)\langle
m_{\alpha}(-{\bf q})m_{\beta}({\bf q})\rangle$ (4)

All of our calculations of neutron scattering responses in [2] employ the
dipole approximation with the form factor being defined by Eq. (3). We
employed analytical approximations for $\langle j_{n}(q)\rangle$ given in
[19].

To go beyond the dipole approximation, requires a more sophisticated treatment
of $\hat{{\bf M}}({\bf q})$. To achieve this, one can represent the matrix
elements of $\hat{{\bf M}}({\bf q})$ between different angular momentum states
$|J,m_{J}\rangle$ as [20, 21]:

$\displaystyle\langle
J,m_{J}|\hat{M}_{x}|J,m_{J}^{\prime}\rangle=-\sqrt{4\pi}\mu_{B}\sum_{K=1,3,5...}\frac{Z(K,q)}{K}P(K,m_{J}^{\prime},m_{J})\big{[}Y_{K-1,m_{J}^{\prime}-m_{J}+1}(\tilde{\bf{q}})\sqrt{(K+m_{J}-m_{J}^{\prime})(K+m_{J}-m_{J}^{\prime}-1)}$
$\displaystyle\qquad\qquad\qquad\qquad-
Y_{K-1,m_{J}^{\prime}-m_{J}-1}(\tilde{\bf{q}})\sqrt{(K-m_{J}+m_{J}^{\prime})(K-m_{J}+m_{J}^{\prime}-1)}\big{]}$
$\displaystyle\langle
J,m_{J}|\hat{M}_{y}|J,m_{J}^{\prime}\rangle=i\sqrt{4\pi}\mu_{B}\sum_{K=1,3,5...}\frac{Z(K,q)}{K}P(K,m_{J}^{\prime},m_{J})\big{[}Y_{K-1,m_{J}^{\prime}-m_{J}+1}(\tilde{\bf{q}})\sqrt{(K+m_{J}-m_{J}^{\prime})(K+m_{J}-m_{J}^{\prime}-1)}$
$\displaystyle\qquad\qquad\qquad\qquad+Y_{K-1,m_{J}^{\prime}-m_{J}-1}(\tilde{\bf{q}})\sqrt{(K-m_{J}+m_{J}^{\prime})(K-m_{J}+m_{J}^{\prime}-1)}\big{]}$
$\displaystyle\langle
J,m_{J}|\hat{M}_{z}|J,m_{J}^{\prime}\rangle=-\sqrt{4\pi}\mu_{B}\sum_{K=1,3,5...}\frac{Z(K,q)}{K}P(K,m_{J}^{\prime},m_{J})Y_{K-1,m_{J}^{\prime}-m_{J}}(\tilde{\bf{q}})\sqrt{(K-m_{J}+m_{J}^{\prime})(K+m_{J}-m_{J}^{\prime})}$
(5)

where $Y_{nl}(\tilde{\bf{q}})$ are spherical harmonic functions of the
direction of ${\bf q}$, $Y_{nl}=0$ for $l>n$ and

$\displaystyle
P(K,m_{J},m_{J}^{\prime})=(-1)^{J-m_{J}^{\prime}}\frac{\begin{pmatrix}K&J&J\\\
m_{J}-m_{J}^{\prime}&m_{J}^{\prime}&-m_{J}\end{pmatrix}}{\begin{pmatrix}K&J&J\\\
0&J&-J\end{pmatrix}}$ (6)

with $\begin{pmatrix}M_{1}&M_{2}&M_{3}\\\ m_{1}&m_{2}&m_{3}\end{pmatrix}$
being Wigner $3j$-symbols. The functions $Z(K,q)$ are linear combinations of
the expectation values of spherical Bessel functions $\langle
j_{n}(q)\rangle$:

$\displaystyle Z(K,q)=c_{K-1}\langle j_{K-1}(q)\rangle+c_{K+1}\langle
j_{K+1}(q)\rangle$ (7)

with the coefficients $c_{K}$ for rare-earth ions being given in Table 6.1 of
[20].

Considering only the $K=1$ term of Eq. (5) results in an equation of the same
form as the dipole approximation [Eq. (2)], with a slightly different form
factor $f^{\prime}(q)$. Using the coefficients $c_{K}$ given for Ce3+ in [20],
we obtain

$\displaystyle\hat{{\bf M}}({\bf q})\approx f^{\prime}(q){\bf m}({\bf q})$ (8)
$\displaystyle f^{\prime}(q)=\langle j_{0}(q)\rangle+\frac{8}{5}\langle
j_{2}(q)\rangle.$ (9)

cf. $f_{\rm dip}(q)$ [Eq. (3)] which differs in that the coefficient in front
of $\langle j_{2}(q)\rangle$ is $\frac{(2-g_{J})}{g_{J}}=\frac{4}{3}$ for
Ce3+.

The terms involving higher values of $K$ contribute higher order multipole
scattering terms, like those discussed in [1]. The $K=3$ and $K=5$ terms come
with associated form factors:

$\displaystyle h(q)=\langle j_{2}(q)\rangle+\frac{10}{3}\langle
j_{4}(q)\rangle$ (10) $\displaystyle g(q)=\langle j_{4}(q)\rangle+12\langle
j_{6}(q)\rangle$ (11)

as found in [1].

Thus, the essential requirement for the validity of the dipole approximation
is that $\langle j_{0}(q)\rangle$ dominates the other functions $\langle
j_{n=2,4,6}(q)\rangle$ over the relevant range of $q$. As long as this is the
case neither the difference between $f^{\prime}(q)$ and $f_{\rm dip}(q)$ nor
the contributions of the $K>1$ terms will be important.

In Fig. 1 we plot the squares of the form factors $f_{\rm
dip}(q),f^{\prime}(q),h(q),g(q)$, since the form factors ultimately appear
quadratically in the scattering intensity. Over the range
$q\in[0,2.5]{\rm\AA}^{-1}$, which is the relevant range for our experiments,
$f_{\rm dip}(q)^{2}$ agrees closely with $f^{\prime}(q)^{2}$ and the higher-
order contributions $h(q)^{2}$ and $g(q)^{2}$ are negligible. This
demonstrates the validity of our use of the dipole approximation in making
theoretical calculations to compare with neutron scattering.

##  III Use of pseudospins in describing physical properties

One of the foremost complaints in [1] is that we used pseudospin-$1/2$
operators to describe the scattering, which the author argues have no physical
relevance. In this section, we will explain the use of these operators and
demonstrate that they are appropriate for the purposes for which we used them.
We note that the use of such operators to describe physical properties of
rare-earth magnets is a well-established method in the literature, with a
history of success, including in the description of scattering experiments [3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13].

As mentioned above, the crystal field interactions in Ce2Zr2O7 separate a
Kramers doublet from the rest of the $J=5/2$ multiplet, by a gap $\Delta_{\rm
CEF}\approx 56{\rm meV}$ [14, 15]. This should be compared to the temperature
scales of our experiment $T\lesssim 10$ K $\implies k_{B}T\lesssim 0.86{\rm
meV}$, the estimated interaction energies $J_{\tilde{\alpha}}<0.07$ meV and to
the maximum energy transfers measured in our inelastic neutron scattering
experiments $\omega\lesssim 0.4\ {\rm meV}$. Given these comparisons, we can
be confident that the higher CEF levels are not significantly thermally
populated in our experiment, interactions do not appreciably mix higher CEF
levels into the ground state and our inelastic experiments do not measure
transitions out of the low energy Kramers doublet. Thus, all of the properties
we measure in the experiment should be understood in terms of operators acting
in the space of states defined by the low energy Kramers doublet. This is
precisely the space in which the pseudospins act.

A basis for the low energy doublet is given by the $|J=5/2,m_{J}=\pm
3/2\rangle$ states of the Ce3+ ions, with the quantization axis
${\bf{\hat{z}}}$ being chosen to correspond with the local axis of $C_{3}$
rotational symmetry. We can define a set of pseudospins in this basis
$S^{\alpha}=\frac{1}{2}\sigma^{\alpha}$, with $\sigma^{\alpha}$ being the
Pauli matrices. When combined with $2\times 2$ identity matrix, these matrices
provide a complete basis of hermitian $2\times 2$ matrices. This means that
all physical processes which involve only the low energy states of the doublet
can be represented as combinations of these pseudospin operators. Due to the
temperature and energy scales of our experiment, this comprises all of the
processes to which we are sensitive. This is why the pseudospin operators are
the appropriate tool for calculations.

Once projected into this low energy basis, all physical properties can be
represented in the pseudospin language. For example, defining $P$ to be the
projector which projects into the low energy space, we see that the magnetic
moment operator is:

$\displaystyle P{\bf m}P=g_{J}\mu_{B}P{\bf J}P=\mu_{B}\begin{pmatrix}0\\\ 0\\\
g_{z}S^{z}\end{pmatrix}$ (12)

where $g_{z}=18/7\approx 2.57$.

Not only the magnetic moment, but higher order multipoles can also be
represented in terms of pseudospins. For example, we consider the octupolar
operator $\overline{J_{x}J_{x}J_{y}}-J_{y}^{3}$ (where $\overline{abc}$
represents symmetrization over operator orderings). Projecting this into the
low energy doublet gives

$\displaystyle P(\overline{J_{x}J_{x}J_{y}}-J_{y}^{3})P=24S^{y}.$ (13)

All other higher-order multipoles, can also be represented in this way once
projected into the low energy doublet. Thus, the pseudospin operators have
simple and direct relations to physically measurable quantities.

We conclude this section of the reply by providing the derivation of Eq. (F1)
of [2] which gives the neutron scattering intensity in terms of the
correlations of pseudospins, within the dipole approximation [Eq. (4)].

The pseudospin basis used in [2] is related to the one defined above, by a
rotation around the $y$-axis of pseudospin space by an angle $\theta$. This
coordinate frame is chosen to remove a coupling between the $x$ and $z$
pseudospin components from the Hamiltonian (see [16, 5]). The operators in the
new coordinate frame $\tilde{x},\tilde{y},\tilde{z}$ are related to the
original one $x,y,z$ as follows:

$\displaystyle
S^{x}_{i}=\cos(\theta)S^{\tilde{x}}_{i}-\sin(\theta)S^{\tilde{z}}_{i}$ (14)
$\displaystyle S^{y}_{i}=S^{\tilde{y}}_{i}$ (15) $\displaystyle
S^{z}_{i}=\cos(\theta)S^{\tilde{z}}_{i}+\sin(\theta)S^{\tilde{x}}_{i}.$ (16)

Combining this with Eq. (12), we obtain the magnetic moment on each site as:

$\displaystyle{\bf
m}_{i}=g_{z}\mu_{B}(\cos(\theta)S^{\tilde{z}}_{i}+\sin(\theta)S^{\tilde{x}}_{i})\hat{{\bf
z}}_{i}$ (17)

where $\hat{{\bf z}}_{i}$ is the local axis of $C_{3}$ symmetry on site $i$.
There are four distinct sites in the unit cell of the pyrochlore lattice,
defining four sublattices, with different anisotropy axes $\hat{{\bf z}}_{i}$.

We define a lattice Fourier transformed operator ${\bf m}_{j}({\bf q})$ as the
Fourier transform over the sublattice $j$, with $j=0,1,2,3$ being a sublattice
index. Then ${\bf m}({\bf q})=\sum_{j}{\bf m}_{j}({\bf q})$. Inserting this
into Eq. (4) we obtain

$\displaystyle S({\bf q})=g_{z}^{2}\mu_{B}^{2}f_{\rm dip}({\bf
q})^{2}\sum_{i,j,=0,1,2,3}\left(\hat{\bf z}_{i}\cdot\hat{\bf
z}_{j}-\frac{(\hat{\bf z}_{i}\cdot{\bf q})(\hat{\bf z}_{j}\cdot{\bf
q})}{q^{2}}\right)\langle(\cos(\theta)S^{\tilde{z}}_{i}(-{\bf
q})+\sin(\theta)S^{\tilde{x}}_{i}(-{\bf
q}))(\cos(\theta)S^{\tilde{z}}_{j}({\bf q})+\sin(\theta)S^{\tilde{x}}_{j}({\bf
q}))\rangle$ (18)

Since we compare only relative, rather than absolute intensities with
experiment we drop the multiplicative constant $g_{z}^{2}\mu_{B}^{2}$. Finally
we note that the symmetries of $XYZ$ Hamiltonian cause the cross correlators
$\langle S^{\tilde{z}}(-{\bf q})S^{\tilde{x}}({\bf q})\rangle$ to vanish.
Applying these points to Eq. (18) we obtain Eq. (F1) of [2]:

$\displaystyle S({\bf q})=f_{\rm dip}({\bf
q})^{2}\sum_{i,j,=0,1,2,3}\left(\hat{\bf z}_{i}\cdot\hat{\bf
z}_{j}-\frac{(\hat{\bf z}_{i}\cdot{\bf q})(\hat{\bf z}_{j}\cdot{\bf
q})}{q^{2}}\right)\left[\cos(\theta)^{2}\langle S^{\tilde{z}}_{i}(-{\bf
q})S^{\tilde{z}}_{j}({\bf q})\rangle+\sin(\theta)^{2}\langle
S^{\tilde{x}}_{i}(-{\bf q})S^{\tilde{x}}_{j}({\bf q}))\rangle\right]$ (19)

Eqs. (E1) and (E2) of [2] follow similarly, with consideration of the neutron
polarization.

##  IV Inconsistency of the comment’s calculations with experiment

Finally, we comment briefly on the calculations of the scattering amplitude in
[1]. In the comment, the expectation value of the scattering operator ${\bf
Q}$ (proportional to $\hat{{\bf M}}$ as given in Eq. (5)) is calculated,
assuming two different ground states.

In both cases the ground state of the many body system is assumed to be
approximated as a product of single site wavefunctions. The two single-site
wavefunctions proposed are:

$\displaystyle|u\rangle=a|J=5/2,m_{J}=3/2\rangle+b|J=5/2,m_{J}=-3/2\rangle$
(20)
$\displaystyle|g\rangle=\frac{1}{\sqrt{2}}\big{(}(a+ib)|J=5/2,m_{J}=3/2\rangle$
$\displaystyle\qquad+(b-ia)|J=5/2,m_{J}=-3/2\rangle\big{)}$ (21)

with $a$ and $b$ being unknown real coefficients and $a^{2}+b^{2}=1$.

It is important to note here that both $|u\rangle$ and $|g\rangle$ break time
reversal symmetry. This can be seen in the case of $|u\rangle$ by the
generally finite expectation value of $J_{z}$

$\displaystyle\langle u|J_{z}|u\rangle=\frac{3}{2}(a^{2}-b^{2})$ (22)

implying magnetic order. As noted in [1] this also breaks some lattice
symmetries. In fact, this would correspond with all-in-all-out
antiferromagnetic order, which we do not observe.

The state $|g\rangle$ preserves all lattice symmetries, but nevertheless
possesses a finite value of the octupolar order parameter

$\displaystyle\langle
g|\left(\overline{J_{x}J_{x}J_{y}}-J_{y}^{3}\right)|g\rangle=12$ (23)

which breaks time reversal symmetry. This corresponds with the octupole
ordered state on the phase diagram.

Since both these ground states break symmetry, reaching them from the high
temperature paramagnet would require a thermodynamic phase transition. No such
phase transition is observed in either our heat capacity or susceptibility
measurements. Muon spin relaxation measurements [14], also find no evidence of
time reversal symmetry breaking. Thus the assumptions behind the calculation
in [1] are inconsistent with experiment.

In the absence of time reversal symmetry breaking, the expectation value
$\langle{\bf Q}\rangle$ vanishes and the correlator $\langle{\bf Q}\cdot{\bf
Q}\rangle$ cannot be approximated as $\langle{\bf Q}\rangle\cdot\langle{\bf
Q}\rangle$ as suggested in [1].

##  V Summary

To summarize, the arguments advanced in [1] do not warrant any revision to our
analysis of the neutron scattering data in [2]. The multipolar corrections
pointed out in [1] are negligible at the wavevectors we study, and the
pseudospin formalism we use is appropriate for the temperature and energy
range of our experiments. Further, the calculations of the scattering matrix
elements in [1] are themselves not consistent with experiment since they
assume time-reversal symmetry breaking order, requiring a thermodynamic phase
transition on cooling from the high temperature paramagnet, which we do not
observe.

## References

  * [1] S. W. Lovesey, Comment on ”case for a u(1)pi quantum spin liquid ground state in the dipole-octupole pyrochlore ce2zr2o7” by e. m. smith et al., phys. rev. x 12, 021015 (2022), arXiv:2209.03235 .
  * Smith  _et al._ [2022] E. M. Smith, O. Benton, D. R. Yahne, B. Placke, R. Schäfer, J. Gaudet, J. Dudemaine, A. Fitterman, J. Beare, A. R. Wildes, S. Bhattacharya, T. DeLazzer, C. R. C. Buhariwalla, N. P. Butch, R. Movshovich, J. D. Garrett, C. A. Marjerrison, J. P. Clancy, E. Kermarrec, G. M. Luke, A. D. Bianchi, K. A. Ross, and B. D. Gaulin, Case for a ${\mathrm{u}(1)}_{\pi}$ quantum spin liquid ground state in the dipole-octupole pyrochlore ${\mathrm{ce}}_{2}{\mathrm{zr}}_{2}{\mathrm{o}}_{7}$, Phys. Rev. X 12, 021015 (2022).
  * Ross  _et al._ [2011] K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, Quantum excitations in quantum spin ice, Phys. Rev. X 1, 021002 (2011).
  * Savary  _et al._ [2012] L. Savary, K. A. Ross, B. D. Gaulin, J. P. C. Ruff, and L. Balents, Order by quantum disorder in ${\mathrm{er}}_{2}{\mathrm{ti}}_{2}{\mathbf{o}}_{7}$, Phys. Rev. Lett. 109, 167201 (2012).
  * Benton [2016] O. Benton, Quantum origins of moment fragmentation in ${\mathrm{nd}}_{2}{\mathrm{zr}}_{2}{\mathrm{o}}_{7}$, Phys. Rev. B 94, 104430 (2016).
  * Li  _et al._ [2015] Y. Li, G. Chen, W. Tong, L. Pi, J. Liu, Z. Yang, X. Wang, and Q. Zhang, Rare-earth triangular lattice spin liquid: A single-crystal study of ${\mathrm{ybmggao}}_{4}$, Phys. Rev. Lett. 115, 167203 (2015).
  * Rau  _et al._ [2016] J. G. Rau, L. S. Wu, A. F. May, L. Poudel, B. Winn, V. O. Garlea, A. Huq, P. Whitfield, A. E. Taylor, M. D. Lumsden, M. J. P. Gingras, and A. D. Christianson, Anisotropic exchange within decoupled tetrahedra in the quantum breathing pyrochlore ${\mathrm{ba}}_{3}{\mathrm{yb}}_{2}{\mathrm{zn}}_{5}{\mathrm{o}}_{11}$, Phys. Rev. Lett. 116, 257204 (2016).
  * Wen  _et al._ [2017] J.-J. Wen, S. M. Koohpayeh, K. A. Ross, B. A. Trump, T. M. McQueen, K. Kimura, S. Nakatsuji, Y. Qiu, D. M. Pajerowski, J. R. D. Copley, and C. L. Broholm, Disordered route to the coulomb quantum spin liquid: Random transverse fields on spin ice in ${\mathrm{pr}}_{2}{\mathrm{zr}}_{2}{\mathrm{o}}_{7}$, Phys. Rev. Lett. 118, 107206 (2017).
  * Rau and Gingras [2019] J. G. Rau and M. J. Gingras, Frustrated quantum rare-earth pyrochlores, Annu. Rev. Condens. Matter Phys 10, 357 (2019).
  * Sarkis  _et al._ [2020] C. L. Sarkis, J. G. Rau, L. D. Sanjeewa, M. Powell, J. Kolis, J. Marbey, S. Hill, J. A. Rodriguez-Rivera, H. S. Nair, D. R. Yahne, S. Säubert, M. J. P. Gingras, and K. A. Ross, Unravelling competing microscopic interactions at a phase boundary: A single-crystal study of the metastable antiferromagnetic pyrochlore ${\mathrm{yb}}_{2}{\mathrm{ge}}_{2}{\mathrm{o}}_{7}$, Phys. Rev. B 102, 134418 (2020).
  * Sibille  _et al._ [2020] R. Sibille, N. Gauthier, E. Lhotel, V. Porée, V. Pomjakushin, R. Ewings, T. Perring, J. Ollivier, A. Wildes, C. Ritter, T. Hansen, D. Keen, G. Nilsen, L. Keller, S. Petit, and T. Fennell, A quantum liquid of magnetic octupoles on the pyrochlore lattice, Nat. Phys. 16, 546–552 (2020).
  * Dun  _et al._ [2020] Z. Dun, X. Bai, J. A. M. Paddison, E. Hollingworth, N. P. Butch, C. D. Cruz, M. B. Stone, T. Hong, F. Demmel, M. Mourigal, and H. Zhou, Quantum versus classical spin fragmentation in dipolar kagome ice ${\mathrm{ho}}_{3}{\mathrm{mg}}_{2}{\mathrm{sb}}_{3}{\mathrm{o}}_{14}$, Phys. Rev. X 10, 031069 (2020).
  * Samartzis  _et al._ [2022] A. Samartzis, J. Xu, V. K. Anand, A. T. M. N. Islam, J. Ollivier, Y. Su, and B. Lake, Pinch points and half-moons in dipolar-octupolar ${\mathrm{nd}}_{2}{\mathrm{hf}}_{2}{\mathrm{o}}_{7}$, Phys. Rev. B 106, L100401 (2022).
  * Gao  _et al._ [2019] B. Gao, T. Chen, D. Tam, C.-L. Huang, K. Sasmal, D. Adroja, F. Ye, H. Cao, G. Sala, M. Stone, C. Baines, J. Verezhak, H. Hu, J.-H. Chung, X. Xu, S.-W. Cheong, M. Nallaiyan, S. Spagna, M. Maple, A. Nevidomskyy, E. Morosan, G. Chen, and P. Dai, Experimental signatures of a three-dimensional quantum spin liquid in effective spin-1/2 Ce2Zr2O7 pyrochlore, Nat. Phys. 15, 1052–1057 (2019).
  * Gaudet  _et al._ [2019] J. Gaudet, E. M. Smith, J. Dudemaine, J. Beare, C. R. C. Buhariwalla, N. P. Butch, M. B. Stone, A. I. Kolesnikov, G. Xu, D. R. Yahne, K. A. Ross, C. A. Marjerrison, J. D. Garrett, G. M. Luke, A. D. Bianchi, and B. D. Gaulin, Quantum spin ice dynamics in the dipole-octupole pyrochlore magnet Ce2Zr2O7, Phys. Rev. Lett. 122, 187201 (2019).
  * Huang  _et al._ [2014] Y.-P. Huang, G. Chen, and M. Hermele, Quantum spin ices and topological phases from dipolar-octupolar doublets on the pyrochlore lattice, Phys. Rev. Lett. 112, 167203 (2014).
  * Benton [2020] O. Benton, Ground-state phase diagram of dipolar-octupolar pyrochlores, Phys. Rev. B 102, 104408 (2020).
  * Patri  _et al._ [2020] A. S. Patri, M. Hosoi, and Y. B. Kim, Distinguishing dipolar and octupolar quantum spin ices using contrasting magnetostriction signatures, Phys. Rev. Research 2, 023253 (2020).
  * Lisher and Forsyth [1971] E. J. Lisher and J. B. Forsyth, Analytic approximations to form factors, Acta Crystallogr., Sect. A. 27, 545 (1971).
  * Marshall and Lovesey [1971] W. Marshall and S. W. Lovesey,  _Theory of Thermal Neutron Scattering_ (Oxford Universiy Press, 1971).
  * Rotter and Boothroyd [2009] M. Rotter and A. T. Boothroyd, Going beyond the dipole approximation to improve the refinement of magnetic structures by neutron diffraction, Phys. Rev. B 79, 140405 (2009).


11institutetext: A*star, Institute of High Performance Computing, Singapore
11email<EMAIL_ADDRESS>22institutetext: York University,
Centre for Vision Research, Canada
22email<EMAIL_ADDRESS>
# A Reliable Online Method for Joint Estimation of Focal Length and Camera
Rotation

Yiming Qian 11 0000-0002-1795-2038 James H. Elder 22

###### Abstract

Linear perspective cues deriving from regularities of the built environment
can be used to recalibrate both intrinsic and extrinsic camera parameters
online, but these estimates can be unreliable due to irregularities in the
scene, uncertainties in line segment estimation and background clutter. Here
we address this challenge through four initiatives. First, we use the
PanoContext panoramic image dataset [27] to curate a novel and realistic
dataset of planar projections over a broad range of scenes, focal lengths and
camera poses. Second, we use this novel dataset and the YorkUrbanDB [4] to
systematically evaluate the linear perspective deviation measures frequently
found in the literature and show that the choice of deviation measure and
likelihood model has a huge impact on reliability. Third, we use these
findings to create a novel system for online camera calibration we call
$f\mathbf{R}$, and show that it outperforms the prior state of the art,
substantially reducing error in estimated camera rotation and focal length.
Our fourth contribution is a novel and efficient approach to estimating
uncertainty that can dramatically improve online reliability for performance-
critical applications by strategically selecting which frames to use for
recalibration.

##  1 Introduction

Online camera calibration is an essential task for applications such as
traffic analytics, mobile robotics, architectural metrology and sports
videography. While intrinsic parameters can be estimated in the lab, these
parameters drift due to mechanical fluctuations. Online estimation of
extrinsic camera parameters is crucial for translating observations made in
the image to quantitative inferences about the 3D scene. If a fixed camera is
employed, extrinsics can potentially be estimated manually at deployment, but
again there will be drift due to mechanical and thermal variations, vibrations
and wind, for outdoor applications. With longer viewing distances, rotational
drift can lead to major errors that may be devastating for performance-
critical tasks, such as judging the 3D distance between a pedestrian and a
car. For PTZ cameras and mobile applications, camera rotation varies over
time. PTZ encoder readings are not always available and become less reliable
over time [22], and for mobile applications IMU data are subject to drift. For
all of these reasons, reliable visual methods for online geometric
recalibration of camera intrinsics and extrinsics are important.

Coughlan and Yuille [3] introduced an approach to online estimation of 3D
camera rotation based on a “Manhattan World” (3-point perspective) assumption,
i.e., that a substantial portion of the linear structure in the image projects
from three mutually orthogonal directions: one vertical and two horizontal.
This assumption can apply quite generally in the built environment and is
relevant for diverse application domains including traffic analytics, mobile
robotics, architecture and sports videography (Fig. 1).

Figure 1:  Example application domains of online camera calibration. We show
here labelling of line segments according to Manhattan direction, as
determined by our novel fR system that simultaneously estimates focal length
and camera orientation.

Subsequent work, reviewed below, has built on this idea to improve accuracy
and incorporate simultaneous estimation of focal length. However, we have
found that these methods vary a lot in accuracy depending upon the scene and
whether focal length is known or estimated. This motivates our attempt to
better understand how these methods generalize to diverse scenes and knowledge
of intrinsics, and to develop methods to estimate uncertainty of individual
estimates.

We address these limitations through four contributions: 1) From the
PanoContext dataset [27] we curate a novel and realistic PanoContext-fR
dataset of planar projections over a broad range of scenes, focal lengths and
camera poses that can be used to evaluate systems for simultaneous estimation
of focal length and camera rotation (tilt/roll). 2) We use this novel dataset
and the YorkUrbanDB [4] to evaluate the linear perspective deviation measures
found in the literature and show that a) the choice of deviation measure has a
huge impact on reliability, and b) it is critical to also employ the correct
likelihood models. 3) We use these findings to create a novel system for
online camera calibration we call $f\mathbf{R}$, and show that it outperforms
the state of the art. 4) We develop a novel approach to estimating uncertainty
in parameter estimates. This is important for two reasons. First, knowledge of
uncertainty can be factored into risk models employed in engineering
applications, co-determining actions designed to mitigate this risk. Second,
especially for cases where parameters are expected to be slowly varying,
systems can sparsely select frames that are more likely to generate trustable
estimates of camera parameters. The source code and PanoContext-fR dataset are
available on GitHub111https://github.com/ElderLab-York-University/OnlinefR.

##  2 Prior work

###  2.1 Image features and deviation measures

To apply linear perspective cues to camera calibration, two key design choices
are 1) what features of the image to use and 2) how to measure the agreement
of these features with hypothesized vanishing points. Coughlan and Yuille’s
original approach [3] employed a likelihood measure on the angular deviation
between the line connecting a pixel to the hypothesized vanishing point (which
we will refer to in the following as the vanishing line and the ray passing
through the pixel in the direction orthogonal to the local luminance gradient.

There have since been diverse attempts to improve on this approach, in large
part by using different features and different measures of agreement.
Deutscher et al. [5] and Košeckà and Zhang [9] first generalized this approach
to allow simultaneous estimation of focal length. While Deutscher et al.
retained the gradient field representation of Coughlan and Yuille, in their
Video Compass system, Košeckà and Zhang switched to using sparser line
segments, and proposed a Gaussian likelihood model based upon the projection
distance from the line passing through a line segment to its associated
vanishing point (Fig. 2 (a)).

|  |  |  |
---|---|---|---|---
(a) | (b) | (c) | (d) | (e)
Figure 2:  Five deviation measures evaluated (see also [25]): (a) distance
from vanishing point to line passing through line segment. (b) angular
deviation between line segment and the vanishing line through the line segment
midpoint. (c) distance from line segment endpoint to vanishing line (d)
distance from line segment endpoint to vanishing line, measured orthogonal to
line segment. (e) angle between the interpretation plane normal and the plane
orthogonal to the vanishing direction [19].

Rother [16] used line segments as well, employing a heuristic version of the
deviation measure similar to Coughlan and Yuille’s (Fig. 2 (b)). Lee et al.
[10] followed a similar approach to simultaneous estimation of focal length
and camera rotation, and have open-sourced their code, allowing us to evaluate
their system below. Denis et al. [4] employed the same angular deviation
measure for edges, but returned to Coughlan and Yuille’s probabilistic
framework, using a learned generalized exponential distribution to model these
deviations.

Tardif [20] continued the practice of using line segments, but proposed an
alternative deviation measure based on the squared distance between segment
endpoints and the vanishing line (Fig. 2 (c)). Wildenauer [21] later adopted
the same measure but, similar to Denis et al. [4], returned to Coughlan and
Yuille’s probabilistic framework, using a Cauchy distribution to model the
signed distance of segment endpoints to the vanishing line.

These three deviation measures were reviewed by Xu et al. [25], who used
intuitive arguments to suggest that each of these measures is flawed. Based
upon their analysis, they proposed a novel deviation measure that also uses
the distance of the segment endpoints to the vanishing line, but measures this
distance orthogonal to the segment, instead of orthogonal to the vanishing
line (Fig. 2 (d)). Their full probabilistic model also favours vanishing
points that lie far from the segment midpoint, along the segment direction.

Hough maps can be used to accurately detect lines in an image. Motivated by
the early work of Collins & Weiss [2], Tal & Elder [19] proposed a
probabilistic Hough method to estimate lines and then estimated camera
rotation based on a probabilistic model of deviation in the Gauss sphere,
specifically, the angular deviation between the interpretation plane normal
and the plane orthogonal to the vanishing direction (Fig. 2(e)).

###  2.2 Benchmarks

The most prevalent dataset used for Manhattan scene analysis is the
YorkUrbanDB dataset [4], which does sample a broad range of scenes but is
limited in that focal length is fixed and roll and tilt vary over a relatively
modest range. Also, the 101 images in the dataset are too few to train a deep
network. More recently there have been efforts [8, 12] to create datasets
large enough to train deep networks by curating random planar projections from
existing panoramic image datasets such as SUN360 [24] and Google StreetView.
This makes it easy to generate enough images to train a deep network, and also
ensures that we have exact ground truth for camera focal length, tilt and roll
angle. Unfortunately, there is no ground truth for camera pan, and so networks
trained on these datasets are unable to recover the full camera rotation ${\bf
R}$.

###  2.3 State-of-the-Art Systems

In this paper we are focused on the problem of online estimation of focal
length $f$ and camera rotation ${\bf R}$. Unfortunately, many of the systems
reviewed above assumed known focal length, and/or did not open-source their
code to allow comparison. The exception, as mentioned, is the line segment
method of Lee et al. [10] (Fig. 2 (b)) which we compare against below.

More recently, Simon et al. [18] use the angular deviation measure in Fig.
2(b) but introduce a preprocessing step to estimate the horizon as a guide to
estimating Manhattan vanishing points. They have open-sourced their code.

Hold-Geoffroy et al. [8] kicked off the use of deep learning for online camera
calibration, using the Sun360 panoramic image dataset [24] to generate planar
projections based on random samplings of focal length, horizon, roll and
camera aspect ratio. They used this dataset to train, validate and test a
network based on a denseNet backbone and three separate heads to estimate the
horizon height and angle and the vertical FOV, from which camera focal length,
roll and tilt can be derived. While the Hold-Geoffroy et al. [8] dataset has
not been released and the method has not been open-sourced, they do provide a
web interface for testing their method on user-provided images, allowing us to
evaluate and compare their approach on public datasets.

Lee et al. [12] followed a similar approach to generating large datasets (see
above) to train a neural geometric parser. Their system takes line segments
detected by LSD [7] as input, first estimating the zenith vanishing point,
selecting segments that are close to vertical, generating vanishing point
hypotheses from pairs of these segments, and then feeding the segments and the
hypotheses into a deep network to score the hypotheses. High-scoring zenith
candidates are then used to generate camera tilt/roll and focal length
hypotheses that are combined with the image and Manhattan line maps as input
to a second network to generate final camera tilt/roll and focal length
estimates. In a more recent paper, Lee et al. [11] have introduced an updated
system called CTRL-C that continues to employ LSD line segments but shifts to
a transformer architecture for estimating the parameters. While the neural
geometric parser has not been open-sourced, CTRL-C has, and we compare with
their approach below.

There are other recent deep learning approaches to estimation of vanishing
points and/or camera rotation [23, 14, 13], but these do not estimate focal
length.

##  3 $f\mathbf{R}$ Method

To develop our $f\mathbf{R}$ method for online camera calibration, we will
assume a camera-centred world frame aligned with the Manhattan structure of
the scene and employ a standard projection model
$\tilde{\mathbf{x}}=\mathbf{KR}\bar{\mathbf{X}}$ where $\tilde{\mathbf{x}}$ is
an image point in homogenous coordinates, ${\mathbf{K}}$ is the camera’s
intrinsic matrix, ${\mathbf{R}}$ is the camera rotation matrix and
$\bar{\mathbf{X}}$ is a 3D world point in augmented coordinates. Unless
otherwise noted, we will assume that any nonlinear distortions in the camera
have been calibrated out in the lab, and that after laboratory calibration the
camera has a central principal point, square pixels and zero skew, leaving a
single intrinsic unknown, the focal length $f$. While these assumptions will
not be met exactly by real cameras in the field, they are the standard
assumptions employed by all of the methods reviewed above and the methods we
compare against below, and in our experience are reasonable approximations for
most higher-quality modern cameras.

The goal of $f\mathbf{R}$ then is to estimate the unknown focal length $f$,
together with the 3D rotation ${\bf R}$ of the camera relative to the
Manhattan frame of the scene. Again, we note that this goes beyond the
capacity of recent deep learning approaches, which are only able to estimate
two of the three rotational degrees of freedom (roll and tilt). Standard
formulae can be used to convert between camera Euler angles (roll, tilt, pan)
and the camera rotation matrix $\mathbf{R}$ [6].

We collect these unknowns into the parameter set $\Psi=\\{f,{\bf R}\\}$. Note
that $\Psi$ completely determines the locations of the three Manhattan
vanishing points: The $i$th vanishing point is given by
$\tilde{\mathbf{x}}_{i}=\mathbf{KR}_{i}$, where ${\mathbf{R}}_{i}$ is the
$i$th column of ${\mathbf{R}}$.

In order to better understand the relative advantages of the various deviation
measures proposed in the literature, and the role of probabilistic modeling,
we will adopt the original mixture model framework first proposed by Coughlan
& Yuille [3] for luminance gradients and later extended to line segments [9],
edges [4] and lines [19]. We will focus on the use of line segments as
features, as there are numerous high-quality open-source line segment
detectors now available [7, 1, 15] and line segments have been shown to be
effective features for Manhattan scene analysis. We will use the MCMLSD line
segment detector [1] for most of the experiments below, but will also evaluate
the advantages and disadvantages of alternative detectors.

###  3.1 Probabilistic model

$f\mathbf{R}$ will assume that each line segment $\vec{l}_{i}$ in the image is
generated by a model $m_{i}\in M$ corresponding to one of four possible
processes: the three Manhattan families of parallel lines or a background
process assumed to generate a uniform distribution of lines. Given camera
parameters $\Psi$, the probability of observing the segment $\vec{l}_{i}$ is
then given by:

$p(\vec{l}_{i}|\Psi)=\sum_{m_{i}\in M}p(\vec{l}_{i}|\Psi,m_{i})p(m_{i})$ (1)

where $p(m_{i})$ is the prior probability that a segment is generated by
process $m_{i}$ and $p(\vec{l}_{i}|\Psi,m_{i})$ is the likelihood of an
observed line generated by $m_{i}$ under camera parameters $\Psi$.

Since, as noted above, the parameters $\Psi$ completely determine the
vanishing points, the likelihoods $p(\vec{l}_{i}|\Psi,m_{i})$ for Manhattan
processes $m_{i}$ can be parameterized by one of the non-negative measures of
deviation between the line $\vec{l}_{i}$ and the corresponding Manhattan
vanishing point (deviation measures a-e reviewed above and summarized in Fig.
2). We will generally use the YorkUrbanDB training partition to fit
exponential models $p(x)=\frac{1}{\lambda}\exp(-x/\lambda)$ for these
deviations, and assume that the background process is uniform on this measure
- see supplemental material for the fits. The one exception is for deviation
measure d, where we employ the central Gaussian model proposed by Xu et al.
[25], which also favours vanishing points that lie far from the segment
midpoint, along the segment direction. Here the dispersion parameter $\sigma$
does not have a clear interpretation in terms of the YorkUrbanDB so we use the
value of $\sigma=1$ pixel recommended by Xu et al. Table 1 summarizes these
parameters.

Table 1: Dispersion parameters for likelihood models. Deviation measure | Model | Param | Horiz | Vert
---|---|---|---|---
a | Exponential | $\lambda$ | 94.46 pix | 17.26 pix
b | Exponential | $\lambda$ | 1.46 deg | 0.57 deg
c | Exponential | $\lambda$ | 0.39 pix | 0.53 pix
d | Gaussian | $\sigma$ | 1.00 pix | 1.00 pix
e | Exponential | $\lambda$ | 0.80 deg | 0.57 deg

The priors $p(m_{i})$ are also estimated from the YorkUrbanDB training set:
45% of lines are expected to arise from the vertical process, 26% from each of
the horizontal processes and 3% from the background process.

Assuming conditional independence over $n$ observed line segments
$\vec{l}_{i}$ and a flat prior over the camera parameters $\Psi$, the optimal
solution should be found by maximizing the sum of log likelihoods, however we
find empirically that we obtain slightly better results if we weight by line
segment length $\left|\vec{l}_{i}\right|$, solving for

$\Psi^{*}=\arg\max_{\Psi}\sum_{i=1,\ldots,n}\left|\vec{l}_{i}\right|\log
p\left(\vec{l}_{i}|\Psi\right)$ (2)

We conjecture that the small improvement achieved by weighting by line length
may derive from a statistical dependence between line length and the
likelihood model: Longer lines may be more accurate.

###  3.2 Parameter search

The $f\mathbf{R}$ objective function (Eqn. 2) is generally non-convex. We opt
for a simple search method that can be used to compare the deviation measures
a-e summarized in Fig. 2 and to assess the role of probabilistic modeling. The
method proceeds in two stages. In Stage 1, we do a coarse $k\times k\times
k\times k$ grid search over the four-dimensional parameter space of $\Psi$,
evaluating Eqn. 2 at each of the $k^{4}$ parameter proposals. Specifically, we
sample uniformly over the range $[-45,+45]$ deg for pan, $[-15,+15]$ deg for
roll, $[-35,+35]$ deg for tilt and $[50,130]$ deg for horizontal FOV. (We also
evaluated a RANSAC approach [21] for Stage 1 but found it be less accurate -
see supplementary material for results.) We then deploy a nonlinear iterative
search (MATLAB fmincon) initialized at each of the the top $l$ of these
$k^{4}$ proposals and constrained to solutions within the ranges above. In the
following we will use $k=8,l=4$.

Note that focal length $f$ and FOV are monotonically related through ${\rm
FOV}=2\arctan\left(\frac{w}{2f}\right)$, where $w$ is the image width. While
we search over FOV we will generally convert to focal length when reporting
results.

###  3.3 Error prediction

One advantage of a probabilistic approach to online camera calibration is the
opportunity to predict when estimates can be trusted. Since we explicitly
model the distribution of deviations between line segments and hypothesized
vanishing points, we might hope to estimate uncertainty in camera parameters
using standard error propagation methods. Unfortunately, we find empirically
that this local method does not lead to very accurate estimates of
uncertainty, perhaps because the derivative of the camera parameters with
respect to the segment orientations at the estimated parameters $\hat{\psi}$
is not very predictive of the derivative at the true camera parameters
$\psi^{*}$ (see supplementary material for details). We therefore propose
instead three easily computable global predictors to inform the level of trust
that should be invested in estimated camera parameters:

1) Number of line segments. We assign each detected segment$\vec{l}_{i}$ to
the most likely generating process $m_{i}$ under our mixture model (Eqn. 1),
count those assigned to each Manhattan direction, and take the minimum, to
reflect the intuition that this will be the weakest link in the global
parameter estimation process. 2) Entropy. This cue takes advantage of the
first stage of our parameter search. We conjecture that more reliable
parameter estimation will be reflected in a more peaked distribution of
likelihoods over the $k^{4}$ parameter proposals evaluated. To capture this
intuition, we normalize the likelihoods into a discrete probability
distribution and compute its entropy. Low entropy distributions are expected
to be more reliable. 3) Likelihood. We compute the mean log likelihood of the
final estimate output from the second nonlinear iterative stage of our search,
normalized by the number of line segments.

To predict camera parameter MAE from these cues we fit a KNN regression model
on the PanoContext-$fR$ training partition - see Section 5.4 for details.

##  4 Datasets

We evaluate $f\mathbf{R}$ and competing methods on two datasets: the
YorkUrbanDB [4] test partition, and our novel PanoContext-fR dataset, curated
from [27].

The YorkUrbanDB test partition consists of 51 indoor and outdoor images.
Camera focal length $f$ was fixed and estimated in the lab. Manhattan line
segments were hand-labelled to allow estimation of the camera rotation $\bf R$
[4]. While the camera was handheld, the variation in camera roll and tilt is
fairly modest: standard deviation of 0.83 deg for roll and 4.96 deg for tilt.

To allow evaluation over a wider range of camera parameters and to assess how
well recent deep networks generalize compared to geometry-driven approaches,
we introduce PanoContext-fR, a curation of planar projections from the
PanoContext dataset [27], which contains 706 indoor panoramic scenes. We
randomly divide this dataset into two equal training and test partitions of
353 scenes each. We will use the test partition to evaluate and compare
deviation measures, probabilistic models and state-of-the-art systems for
camera parameter estimation. While our fR system is not a machine learning
method and does not require the training partition, we will explore methods
for learning to predict estimation error in Section 3.3 and will use the
training dataset there.

For each of these scenes we generated 15 planar projections, using a standard
$640\times 480$ pixel resolution. Three images were sampled for each of 5
fixed horizontal FOVs from 60 to 120 deg in steps of 15 deg. Pan, roll and
tilt were randomly and uniformly sampled over [-180,+180], [-10,+10] and
[-30,+30] ranges respectively (Fig. 3). This generated 5,295 images for each
of the training and test partitions. This dataset is not intended for network
training, but rather for evaluating generalization of pre-trained and
geometry-driven models.

We sampled FOVs discretely for PanoContext-fR so we could clearly see whether
algorithms are able to estimate focal length (see below). However, we have
also created and will make publicly available an alternative PanoContext-ufR
that samples horizontal FOV randomly and uniformly over [60,120] deg. (See
supplementary material.) Since the the YorkUrbanDB provides the ground truth
rotation matrix, we can evaluate the frame angle error, i.e., the magnitude of
the rotation required to align the estimated camera frame with the ground
truth frame. Since for PanoContext-fR we do not have ground truth pan, we
report mean absolute roll and tilt errors.

Figure 3: Three example planar projections from our new PanoContext-fR
dataset.

##  5 Experiments

###  5.1 Evaluating deviation measures

We begin by applying the $f\mathbf{R}$ search method (Section 3.2) to identify
the parameters $\Psi$ maximizing agreement between line segments and vanishing
points, based on the five deviation measures a-e summarized in Fig. 2. Without
probabilistic modeling, we employed the objective function
$\sum_{i}\log\delta_{i}$, where $\delta_{i}=d_{i},\theta_{i},\Delta\theta_{i}$
depending on which of the five deviation measures is employed. (We found
empirically that logging the deviations before summing improved results.). We
found that none of the deviation measures yielded good results: Mean focal
length error remained above 11% on the York Urban DB and above 34% on the
PanoContext-fR dataset (See Supplementary Material for detailed results.)

Probabilistic modeling of the deviation measure greatly improves results (see
Table 2 for results on the PanoContext-fR dataset and supplementary material
for results on the YorkUrbanDB). We note that this represents the learning of
only 4 parameters (horizontal and vertical dispersions and priors). Moreover
this learning seems to generalize well, as the parameters used to evaluate on
the PanoContext-fR dataset were learned on the YorkUrbanDB training set.

Table 2: Evaluation of deviation measures on the PanoContext-fR training set. Numbers are mean$\pm$standard error.  Dev. measure | Roll MAE (deg) | Tilt MAE (deg) | Focal length MAE (%)
---|---|---|---
a | 4.53$\pm$0.045 | 14.2$\pm$0.14 | 35.6$\pm$0.35
b ($f\mathbf{R}$) | 0.78$\pm$0.02 | 1.59$\pm$0.06 | 8.4$\pm$0.26
c | 0.90$\pm$0.009 | 2.26$\pm$0.022 | 14.1$\pm$0.14
d | 0.81$\pm$0.008 | 1.67$\pm$0.015 | 10.1$\pm$0.08
e | 0.89$\pm$0.009 | 2.13$\pm$0.021 | 19.5$\pm$0.19

Our second observation is that performance also depends strongly on the
deviation measure, even when carefully modeling the likelihood. The large
errors produced by measure a likely reflect the fact that vanishing points are
often far from the principal point, resulting generally in very unpredictable
deviations from the line passing through an associated line segment. The other
measures yield better performance, but we note that deviation measure e tends
to have higher errors for focal length. We believe this is because the Gauss
Sphere measure of deviation suffers from a degeneracy: As focal length tends
to infinity, the interpretation plane normals collapse to a great circle
parallel to the image plane. Thus a hypothesized vanishing point near the
principal point will always generate small deviations, leading to a large
likelihood.

This leaves deviation measures b-d. Method d performs less well on the
YorkUrbanDB dataset and method c has higher errors on the PanoContext-fR
dataset. Overall, method b, the method originally used by Coughlan & Yuille
for isophotes [3], Rother for line segments [16] and Denis et al. for edges
[4], appears to generate the most consistent performance, and so we adopt this
method as our standard $f\mathbf{R}$ deviation measure for the remainder of
the paper.

###  5.2 Evaluating line segment detectors

All of the experiments above employed the MCMLSD line segment detector [1].
Table 3 assesses the sensitivity of our system to this choice by substituting
two alternative detectors: The well-known LSD detector [7] and a more recent
deep learning detector called HT-LCNN [15]. Performance on the PanoContext-fR
dataset is significantly better using MCMLSD than the more recent HT-LCNN,
despite the fact that HT-LCNN is reported [15] to have much better precision-
recall performance on the YorkUrban DB and Wireframe [26] datasets.

We believe this discrepancy stems largely from the incompleteness of these
datasets, which do not claim to label all Manhattan line segments. Methods
like MCMLSD, which attempt to detect all line segments (Manhattan and non-
Manhattan), thus tend to achieve lower precision scores on these datasets,
while deep learning methods like HT-LCNN learn to detect only the labelled
segments and thus achieve higher precision scores. Our experiment reveals that
this higher precision does not translate into better performance but rather
worse performance, presumably because HT-LCNN fails to generate Manhattan
segments that may be unlabelled but are nevertheless useful for estimating the
camera parameters. We proceed with the MCMLSD line segment detector but will
return to this choice when considering run-time efficiency (Section 5.5).

Table 3: Evaluating the choice of line segment detector on the PanoContext-fR training set. Numbers are mean $\pm$ standard error. Methods | Run time | Roll MAE | Tilt MAE | Focal length MAE
---|---|---|---|---
| (sec) | (deg) | (deg) | (%)
MCMLSD [1] | 4.23 | ${\bf 0.78\pm 0.02}$ | ${\bf 1.59\pm 0.06}$ | $\bf{8.4\pm 0.26}$
LSD [7] | 0.40 | $1.23\pm 0.012$ | $3.35\pm 0.03$ | $8.6\pm 0.11$
HT-LCNN [15] | 2.91 | $0.93\pm 0.009$ | $1.91\pm 0.02$ | 10.1$\pm$0.14

###  5.3 Comparison with state of the art

We compare our fR system against the four state-of-the-art systems [10, 18, 8,
11] reviewed in Sec. 2. Three of these [10, 18, 11] are open-sourced and we
downloaded the code. The fourth (Hold-Geoffroy) [8] provides a web interface.
We test two versions of the CTRL-C [11] network: CTRL-C S360 was trained on a
dataset curated from the SUN360 dataset [24]. CTRL-C GSV was trained on a
dataset curated from the Google Street View dataset [12].

Table 4)(top) shows results of this comparison on the York UrbanDB test set.
We report roll and tilt error instead of frame error, since the deep learning
methods (Hold-Geoffroy [8] and CTRL-C [11]) do not estimate pan angle. (See
supplemental material for comparison of pan angle error with the methods of
Lee et al. [10] and Simon et al. [18]). (The Simon et al. method returned a
valid result for only 69% of the York UrbanDB test images and 46% of the
PanoContext-$f\mathbf{R}$ test images; we therefore report their average
performance only for these.)

Table 4: Performance comparison with SOA on the York UrbanDB (top) and
PanoContext-fR (bottom) test sets. Numbers are mean $\pm$ standard error.

YorkUrbanDB | Run time | Roll | Tilt | FOV | Focal len.
---|---|---|---|---|---
| (sec) | (MAE deg) | MAE (deg) | MAE (%) | MAE (%)
Lee[10] | 1.13 | 0.80$\pm$0.2 | 1.33$\pm$0.2 | $9.2\pm 1.15$ | $11.5\pm 1.62$
Simon[18] | 0.44 | 2.63$\pm$0.6 | 11.0$\pm$1 | $11.4\pm 2.56$ | $15.1\pm 2.08$
Hold-Geoffroy[8] | n/a | 1.47$\pm$0.2 | 7.57$\pm$0.6 | $48.7\pm 2.11$ | $37.7\pm 7.02$
CTRL-C[11] S360 | 0.32 | 1.38$\pm$0.2 | 1.09$\pm$0.1 | $29.0\pm 2.50$ | $24.7\pm 1.80$
CTRL-C [11] GSV | 0.32 | 1.37$\pm$0.2 | 3.90$\pm$0.3 | $69.4\pm 1.66$ | $48.8\pm 0.77$
fR | 6.40 | 0.50$\pm$0.1 | 1.16$\pm$0.1 | $\bf{3.8\pm 0.61}$ | $\bf{4.6\pm 0.78}$
Fast-fR | 0.89 | 1.13$\pm$0.4 | 1.62$\pm$0.3 | $5.1\pm 1.04$ | $5.7\pm 1.21$

PanoContext-fR | Roll | Tilt | FOV | Focal len.
---|---|---|---|---
| MAE (deg) | MAE (deg) | MAE (%) | (MAE %)
Lee[10] | 2.02$\pm$ 0.02 | 3.35$\pm$0.03 | $13.4\pm 0.22$ | $22.3\pm 0.22$
Simon[18] | 1.50$\pm$0.02 | 3.18$\pm$ 0.05 | $32.8\pm 0.46$ | $54.0\pm 0.53$
Hold-Geoffroy[8] | 1.58$\pm$0.02 | 3.69$\pm$0.04 | $16.1\pm 0.12$ | $20.5\pm 0.20$
CTRL-C[11] S360 | 1.03$\pm$0.01 | 2.19$\pm$0.02 | $8.1\pm 0.07$ | $13.1\pm 0.13$
CTRL-C [11] GSV | 2.24$\pm$0.02 | 8.98$\pm$ 0.09 | $17.3\pm 0.11$ | $21.9\pm 0.31$
fR | 0.78$\pm$0.02 | 1.59$\pm$0.06 | 5.3$\pm$0.15 | 8.4$\pm$0.26
Fast-fR | 0.89$\pm$0.02 | 1.90$\pm$0.07 | $5.4\pm 0.16$ | $\bf{7.6\pm 0.23}$

All methods have fairly low error in roll, but remember that in the York
UrbanDB, variation in camera roll is limited (standard deviation of 0.83 deg).
Our fR system improves on the next best method (Lee et al. [10]) by 38%. For
tilt estimation, variability in performance across systems is more pronounced,
but three systems (our fR system, Lee et al. [10] and CTRL-C S360 [11] all
perform well. Focal length estimation is the big differentiator. Our fR system
performs much better than all of the other systems, beating the next best
system (Lee et al. [10]) by 57% (reducing MAE from 11.5% to 4.9%). Notice that
the deep learning methods (Hold-Geoffroy [8] and CTRL-C [11]) do not perform
well.

Table 4(bottom) compares these systems on our novel PanoContext-fR dataset.
(While FOV and focal length are monotonically related, we show MAE for both
for the reader’s convenience.) We find that our fR system performs
significantly better that all of the other systems on all three benchmarks
(roll, tilt and focal length), improving on the next best system (CTRL-C S360
[11]) by 23%, 25% and 33%, respectively. (The relatively strong performance of
the CTRL-C S360 system here may derive from the fact that the PanoContext
dataset was original drawn from the SUN360 dataset on which CTRL-C S360 is
trained.)

Fig. 4 compares distributions of estimated and ground truth parameters for the
PanoContext-fR dataset. The distributions for camera rotation are generally
reasonable, although the method of Simon et al. [18] exhibits odd preferences
for specific roll and tilt angles, the Hold-Geoffroy system [8] has an overly
strong bias to zero roll, and when trained on GSV, the CTRL-C system [11]
develops a bias against small upward tilts.

The sampling of five discrete FOVs in the PanoContext-fR dataset allows us to
visualize the sensitivity of the algorithms to focal length. While the
geometry-based algorithms all show modulation with the sampled FOVs (albeit
weak for Simon et al [18]), the deep learning methods do not, each forming a
single broad peak. Our fR system appears to be best able to pick up the FOV
signal from the data without a strong prior bias.

The ability of the deep learning methods to estimate FOV/focal length might in
part be cue to limitations in the range of FOVs in their training datasets.
However, in the supplementary material we show that this does not fully
account for the superiority of our fR system, which we believe derives from
more direct geometry and probabilistic modeling.

|  |  |
---|---|---|---
Figure 4: Distribution of ground truth and estimated camera parameters for the
PanoContext-fR test sets.

###  5.4 Predicting Reliability

We employ three global cues to predict camera parameter estimation error: 1)
The minimum number of segments over the three Manhattan directions, 2) Entropy
over our parameter grid search and 3) mean log likelihood of the final
parameter estimate: See supplementary material for a visualization of how
parameter error varies as a function of these cues. To estimate the
reliability of fR estimates at inference, we use the PanoContext-fR training
set to fit a regression model that uses these three cues jointly to predict
the absolute error in each of the camera parameters: We collect these cues
into a 3-vector, use the training data to compute a whitening transform and
then use KNN regression in the 3D whitened space to estimate error, again
selecting K by 5-fold cross-validation.

Table 5 assesses the ability of this model to predict error on the held-out
PanoContext-fR test set. The table shows the mean error obtained if we accept
only the top 25%, 50%, 75% or 100% of estimates predicted to have least error.
The model is remarkably effective. For example, by accepting only the 25% of
estimates predicted to be most reliable, mean error declines by 24% for roll,
52% for tilt and 66% for focal length. This ability to predict estimation
error can be extremely useful in applications where the camera can be
recalibrated on sparse frames, or when conservative actions can be taken to
mitigate risks.

Table 5: Parameter MAE for 25%, 50%, 75% and 100% of held-out PanoContext-fR test images predicted to be most reliable. % of images | Roll MAE (deg) | Tilt MAE (deg) | Focal length MAE (%)
---|---|---|---
25% | 0.56 | 0.68 | 2.8%
50% | 0.59 | 0.75 | 3.9%
75% | 0.63 | 0.86 | 5.2%
100% | 0.74 | 1.42 | 8.3%

###  5.5 Run time

We ran our experiments on a 2.3GHz 8-Core Intel i9 CPU with 32 GB RAM and an
NVIDIA Tesla V100-32GB GPU. Our fR system as implemented is slower than
competing systems (Table 4), primarily due to our line segment detection
system MCMLSD, which consumes on average 4.23 sec per image (Table 3).

To address this issue, we have created a more efficient version of our system
we call Fast-fR, through two innovations. First, we replace MCMLSD with LSD,
which consumes only 40 msec per image on average. Second, we limit the number
of iterations in our second search stage to 10. This reduces average run time
from 6.4 sec to 0.83 sec with only a modest decline in accuracy (Table 4).

##  6 Limitations

fR relies on the Manhattan World assumption, i.e., that images contain aligned
rectilinear structure. The superior performance of fR on the York Urban DB and
PanoContext-fR datasets, which are representative of common outdoor and indoor
built environments, attests to the real-world applicability of this
assumption. Nevertheless, the method will typically fail on scenes that do not
conform to the Manhattan assumption. This may include images with insufficient
structure (see Supplementary Material for examples), nature scenes and non-
Manhattan built environments (e.g., the Guggenheim museums in New York or
Bilbao). Generalizing the fR method to exploit additional scene regularities,
e.g., Atlanta World [17] or quadrics, is an interesting direction for future
research.

Like all SOA methods reviewed and evaluated here, our fR system assumes no
non-linear distortions, a central principal point, square pixels and zero
skew. These are reasonable approximations, since deviations from these
assumptions can be calibrated out in the factory or lab, and subsequent drift
of these parameters is typically minor relative to changes in focal length and
rotation.

While the Manhattan constraint can also be used to estimate the principal
point, we find that more accurate estimates of focal length and rotation are
obtained by assuming a central principal point. Simultaneous online estimation
of focal length, principal point and camera rotation, as well as radial
distortion, thus remains an interesting and challenging direction for future
research.

##  7 Conclusions

We have employed existing and novel datasets to assess the reliability of
online systems for joint estimation of focal length and camera rotation. We
show that the reliability of geometry-driven systems depends profoundly on the
deviation measure employed and accurate probabilistic modeling. Based on these
findings, we have proposed a novel probabilistic geometry-driven approach
called fR that outperforms four state-of-the-art competitors, including two
geometry-based systems and two deep learning systems. The fR advantage is most
pronounced for estimation of focal length, where the deep learning systems
seem to struggle. We note that the deep learning systems also do not predict
pan angle, and so do not fully solve for the camera rotation. We further
demonstrate an ability to estimate the reliability of specific fR predictions,
leading to substantial gains in accuracy when systems can choose to reject
parameter estimates judged to be unreliable. We believe this will be useful in
many real-world applications, particularly in mobile robotics and autonomous
vehicles.

##  8 Acknowledgement

This work was supported by the University of Toronto, NSERC, the Ontario
Research Fund, the York University VISTA and Research Chair programs (Canada),
and the Agency for Science, Technology and Research (A*STAR) under its RIE2020
Health and Biomedical Sciences (HBMS) Industry Alignment Fund Pre-Positioning
(IAF-PP) Grant No. H20c6a0031 and AI3 HTPO Seed Fund (C211118014) (Singapore).

## References

  * [1] Almazan, E.J., Tal, R., Qian, Y., Elder, J.H.: MCMLSD: A dynamic programming approach to line segment detection. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 2031–2039 (2017)
  * [2] Collins, R.T., Weiss, R.S.: Vanishing point calculation as a statistical inference on the unit sphere. In: International Conference on Computer Vision. vol. 90, pp. 400–403. Citeseer (1990)
  * [3] Coughlan, J.M., Yuille, A.L.: Manhattan world: Compass direction from a single image by Bayesian inference. In: International Conference on Computer Vision. vol. 2, pp. 941–947. IEEE (1999)
  * [4] Denis, P., Elder, J.H., Estrada, F.J.: Efficient edge-based methods for estimating manhattan frames in urban imagery. In: European Conference on Computer Vision. pp. 197–210. Springer (2008)
  * [5] Deutscher, J., Isard, M., MacCormick, J.: Automatic camera calibration from a single Manhattan image. In: European Conference on Computer Vision. pp. 175–188. Springer (2002)
  * [6] Faugeras, O.: Three-Dimensional Computer Vision: A Geometric Viewpoint. MIT Press, Cambridge, MA (1993)
  * [7] von Gioi, R.G., Jakubowicz, J., Morel, J.M., Randall, G.: Lsd: A fast line segment detector with a false detection control. IEEE Transactions on Pattern Analysis & Machine Intelligence (4), 722–732 (2008)
  * [8] Hold-Geoffroy, Y., Sunkavalli, K., Eisenmann, J., Fisher, M., Gambaretto, E., Hadap, S., Lalonde, J.F.: A perceptual measure for deep single image camera calibration. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 2354–2363 (2018)
  * [9] Košecká, J., Zhang, W.: Video compass. In: Heyden, A., Sparr, G., Nielsen, M., Johansen, P. (eds.) European Conference on Computer Vision. pp. 476–490. Springer Berlin Heidelberg, Berlin, Heidelberg (2002)
  * [10] Lee, D.C., Hebert, M., Kanade, T.: Geometric reasoning for single image structure recovery. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 2136–2143. IEEE (2009)
  * [11] Lee, J., Go, H., Lee, H., Cho, S., Sung, M., Kim, J.: CTRL-C: Camera calibration transformer with line-classification. In: International Conference on Computer Vision. pp. 16228–16237 (2021)
  * [12] Lee, J., Sung, M., Lee, H., Kim, J.: Neural geometric parser for single image camera calibration. In: European Conference on Computer Vision. pp. 541–557. Springer (2020)
  * [13] Li, H., Chen, K., Kim, P., Yoon, K.J., Liu, Z., Joo, K., Liu, Y.H.: Learning icosahedral spherical probability map based on bingham mixture model for vanishing point estimation. In: International Conference on Computer Vision. pp. 5661–5670 (2021)
  * [14] Li, H., Zhao, J., Bazin, J.C., Chen, W., Liu, Z., Liu, Y.H.: Quasi-globally optimal and efficient vanishing point estimation in Manhattan world. In: International Conference on Computer Vision. pp. 1646–1654 (2019)
  * [15] Lin, Y., Pintea, S.L., van Gemert, J.C.: Deep hough-transform line priors. In: European Conference on Computer Vision. pp. 323–340. Springer (2020)
  * [16] Rother, C.: A new approach to vanishing point detection in architectural environments. Image and Vision Computing 20(9-10), 647–655 (2002)
  * [17] Schindler, G., Dellaert, F.: Atlanta world: An expectation maximization framework for simultaneous low-level edge grouping and camera calibration in complex man-made environments. In: Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004\. vol. 1, pp. I–I. IEEE (2004)
  * [18] Simon, G., Fond, A., Berger, M.O.: A Simple and Effective Method to Detect Orthogonal Vanishing Points in Uncalibrated Images of Man-Made Environments. In: Bashford-Rogers, T., Santos, L.P. (eds.) EuroGraphics - Short Papers. The Eurographics Association (2016). https://doi.org/10.2312/egsh.20161008
  * [19] Tal, R., Elder, J.H.: An accurate method for line detection and Manhattan frame estimation. In: Park, J.I., Kim, J. (eds.) Asian Conference on Computer Vision-Workshops. pp. 580–593. Springer Berlin Heidelberg, Berlin, Heidelberg (2012)
  * [20] Tardif, J.P.: Non-iterative approach for fast and accurate vanishing point detection. In: IEEE International Conference on Computer Vision. pp. 1250–1257. IEEE (2009)
  * [21] Wildenauer, H., Hanbury, A.: Robust camera self-calibration from monocular images of Manhattan worlds. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 2831–2838. IEEE (2012)
  * [22] Wu, Z., Radke, R.: Keeping a pan-tilt-zoom camera calibrated. Pattern Analysis and Machine Intelligence, IEEE Transactions on 35(8), 1994–2007 (Aug 2013). https://doi.org/10.1109/TPAMI.2012.250
  * [23] Xian, W., Li, Z., Fisher, M., Eisenmann, J., Shechtman, E., Snavely, N.: UprightNet: geometry-aware camera orientation estimation from single images. In: International Conference on Computer Vision. pp. 9974–9983 (2019)
  * [24] Xiao, J., Ehinger, K.A., Oliva, A., Torralba, A.: Recognizing scene viewpoint using panoramic place representation. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 2695–2702. IEEE (2012)
  * [25] Xu, Y., Oh, S., Hoogs, A.: A minimum error vanishing point detection approach for uncalibrated monocular images of man-made environments. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 1376–1383 (2013)
  * [26] Xu, Z., Shin, B., Klette, R.: A statistical method for line segment detection. Computer Vision and Image Understanding 138, 61–73 (2015)
  * [27] Zhang, Y., Song, S., Tan, P., Xiao, J.: PanoContext: A whole-room 3D context model for panoramic scene understanding pp. 668–686 (2014)