equation(#1)

# On the Role of Entropy-based Loss for Learning Causal Structures with
Continuous Optimization

Ruichu Cai 1, 2, Weilin Chen1, Jie Qiao1, 2 , Zhifeng Hao1, 3
1 School of Computer Science, Guangdong University of Technology, China
2 Pazhou Lab, Guangzhou, China.
3 School of Mathematics and Big Data, Foshan University, China
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Corresponding author.

###### Abstract

Causal discovery from observational data is an important but challenging task
in many scientific fields. Recently, NOTEARS (Zheng et al., 2018) formulates
the causal structure learning problem as a continuous optimization problem
using least-square loss with an acyclicity constraint. Though the least-square
loss function is well justified under the standard Gaussian noise assumption,
it is limited if the assumption does not hold. In this work, we theoretically
show that the violation of the Gaussian noise assumption will hinder the
causal direction identification, making the causal orientation fully
determined by the causal strength as well as the variances of noises in the
linear case and the noises of strong non-Gaussianity in the nonlinear case.
Consequently, we propose a more general entropy-based loss that is
theoretically consistent with the likelihood score under any noise
distribution. We run extensive empirical evaluations on both synthetic data
and real-world data to validate the effectiveness of the proposed method and
show that our method achieves the best in Structure Hamming Distance, False
Discovery Rate, and True Positive Rate matrices.

##  1 Introduction

Learning causal structure from observational data has become an important
topic in many scientific fields, such as economics (Ghysels et al., 2016),
biology (Grosse-Wentrup et al., 2016), and social science (Cai et al., 2016).
Due to the expensive cost or the ethic of randomized experiments, the task of
causal discovery from purely observational data has been drawn much attention.

Many approaches have been proposed for learning causal structure.
Traditionally, by utilizing the conditional independence property among
observed variables, the constraint-based approaches have been proposed, e.g.,
PC algorithm (Spirtes et al., 2000), but only identify the underlying directed
acyclic graph (DAG) up to Markov equivalence class (Andersson et al., 1997).
Alternatively, by introducing a certain class of Structure Causal Model (SCM),
and further assuming the causal mechanism that the noise and the hypothetical
cause are independent, the functional-based causal models have been proposed,
e.g., the Linear Non-Gaussian Acyclic Model (LiNGAM) (Shimizu et al., 2006),
the Additive Noise Model (ANM) (Hoyer et al., 2008), and the Post-Nonlinear
causal Model (PNL) (Zhang and Hyvärinen, 2009). However, due to the
intractable search space superexponential in graph nodes, learning DAGs using
functional-based causal models is challenging.

Recently, NOTEARS (Zheng et al., 2018), under the additive noise model
assumption, formulates the causal structure learning problem as a continuous
optimization problem using least-square loss with a continuous DAG constraint.
Such a technique has been extensively developed and applied to learning linear
or nonlinear causal structures. Yu et al. (2019) introduces a variational
autoencoder framework for modeling the generative process of a causal
structure equipped with evidence lower bound with a Gaussian prior of noise
which is implemented by least-square loss. Ng et al. (2019b) and Zheng et al.
(2020) extend a linear causal model into a nonlinear causal model using neural
networks but still rely on the least-square loss for reconstruction.

However, our analysis shows that using least-square loss as the score function
will hinder the causal direction identification. In a linear system, such loss
will be highly influenced by the causal strengths and the noise variances
leading to a bias estimation comparing with the likelihood score. In a more
general nonlinear system, the causal direction will be incorrectly identified
if noise distribution has strong non-Gaussianity property.

In this paper, we show that the entropy-based loss is consistent with the
likelihood score under the additive noise model using any noise distribution,
and thus we advocate using the entropy-based loss instead of least-square
loss. Overall, our contributions are as follows:

  1. 1.

In Section 4, we study the limitations of least-square loss and provide
theoretical bounds that least-square loss will fail to identify the causal
direction in a linear and nonlinear system, respectively.

  2. 2.

In Section 5, we build a connection between the entropy-based loss and the
log-likelihood score, and further provide the theory for its validness.

  3. 3.

In Section 6, we run extensive experiments over linear and nonlinear systems
by using both synthetic data and real-world data and show our method performs
more effectively and stably.

##  2 Related Works

Related to least-square loss: Blöbaum et al. (2018, 2019) address the problem
of inferring the causal relation between two variables using least-square loss
under the modularity property, i.e., the independence between the function and
the distribution of cause. Peters and Bühlmann (2014) prove the full
identifiability of the linear Gaussian structural equation in case that all
noise variables have the same variance.

Related to NOTEARS: Yu et al. (2019) proposes an alternative characterization
of acyclicity and utilizes a generative model to learn the nonlinear causal
structure. Ng et al. (2019b) utilizes a graph autoencoder framework to extend
the linear SCM to nonlinear SCM. Ng et al. (2019a) extends to nonlinearity by
putting the weighted matrix on the first layer of MLP. Lachapelle et al.
(2019) expresses weighted matrix by neural network paths to deal with the
nonlinear case. Zheng et al. (2020) proposes a more general acyclicity
constraint based on partial derivatives to support nonlinear models.

##  3 Preliminary

###  3.1 Structure Causal Model

Causal discovery problems can be formalized by a Structure Causal Model (SCM)
(Didelez and Pigeot, 2001). Given a set of random variables
$X=\left\\{X_{1},X_{2},\dots,X_{d}\right\\}$ and the corresponding noises
$N=\left\\{N_{1},N_{2},\dots,N_{d}\right\\}$, SCM between $X_{i}$ and its
direct parents $X_{pa_{\mathcal{G}}\left(i\right)}$ with respect to a DAG
$\mathcal{G}$ is defined as follows:

$X_{i}=f_{i}\left(X_{pa_{\mathcal{G}}\left(i\right)},N_{i}\right),$ (1)

where $f_{i}$ could be a linear or nonlinear function. Following the previous
works (Zheng et al., 2018), we consider the additive noise model, i.e.,
$X_{i}=f_{i}\left(X_{pa_{\mathcal{G}}\left(i\right)}\right)+N_{i}$, and assume
the causal sufficiency holds.

###  3.2 NOTEARS

Zheng et al. (2018) proposed a continuous optimization for learning causal
structure, using least-square loss with an acyalicity constraint. Given a set
of observational data
$\left\\{x_{1}^{(i)},\dots,x_{d}^{(i)}\right\\}_{i=1}^{m}$ which is sampled
from a joint distribution $p(X)$, we assume that the joint ditritbution is
Markov with respect to a ground truth DAG $\mathcal{G}$ and can be factorized
as $p(X)=\prod_{i=1}^{d}p(X_{i}|X_{pa_{\mathcal{G}}(i)})$. Then the
corresponding design matrix is $\mathbf{X}\in\mathbb{R}^{m\times d}$. In
particular, the directed graph in a linear SCM can be encoded by a weighted
adjacency matrix $W\in\mathbb{R}^{d\times d}$, i.e., $X=W^{\mathsf{T}}X+N$.
Zheng et al. (2018) shows that $W$ represents a DAG if and only if
$tr\left(e^{W\circ W}\right)-d=0$ holds, where $\circ$ denotes the Hadamard
product. Then, NOTEARS formulates the causal structure leaning problem as the
following continuous optimization problem:

$\mathop{\arg\min}_{W\in\mathbb{R}^{d\times
d}}\frac{1}{2m}\|\mathbf{X}-\mathbf{X}W\|^{2}_{F}+\lambda\left\|W\right\|_{1}\quad\text{subject
to}\quad tr\left(e^{W\circ W}\right)-d=0,$ (2)

where $\frac{1}{2m}\|\mathbf{X}-\mathbf{X}W\|^{2}_{F}$ is the least-square
loss and is equal, up to constant, to the log-likelihood score of a linear
Gaussian DAG with equal noise variances, and $\left\|W\right\|_{1}$ denotes
the $\ell_{1}$ penalty term on the causal structure.

##  4 Limitations of Least-Square Loss

(a) Accuracy varies with various causal strength

(b) Accuracy varies with various variance of $N_{X}$

(c) Accuracy varies with various variance of $N_{Y}$

Figure 1: Control experiments of accuracy with respect to different conditions
in linear system. At each experiments, we will control one of the parameters
while fix others. The default parameters are set as $\alpha=0.5$,
$\sigma_{N_{X}}=2$, and $\sigma_{N_{Y}}=1$. And the controlled parameters will
varies form $\alpha=[0.5,1]$, $\sigma_{N_{X}}=[0.6,1.6]$,
$\sigma_{N_{Y}}=[0.6,1.6]$.

This section discusses why methods based on least-square loss fail to discover
the underlying causal structure, and we theoretically show the conditions that
they fail in both linear and nonlinear cases.

Without loss of generality, we analyze the causal pair of two variables, and
assume that:

$\displaystyle X$ $\displaystyle=N_{X},$ (3) $\displaystyle Y$
$\displaystyle=\begin{cases}\alpha X+N_{Y}&\text{Linear}\\\
f(X)+N_{Y}&\text{Nonlinear}\end{cases}$

where the variance $\text{Var}(N_{X})=\sigma_{N_{X}}^{2}$,
$\text{Var}(N_{Y})=\sigma_{N_{Y}}^{2}$, $N_{X}\Vbar N_{Y}$, and $f$ is a
nonlinear function.

To see the limitations of least-square loss, we take the two variables as an
example. We run NOTEARS 100 times using the least-square loss with 400 samples
per experiment in linear Gaussian, Uniform, and Gumbel data. As shown in
Figure 1, the accuracy is highly controlled by the causal strength as well as
the noise variance. We see that using the least-square loss, NOTEARS could
fail even under the non-Gaussian distribution that is notably identifiable
(Shimizu et al., 2006). The reason is the mismatch between the truth
distribution and the standard Gaussian distribution implied by the least-
square loss. In addition, due to the non-identifiability of the linear
Gaussian causal model, the causal direction can still be controlled by the
causal strength and the noise variance.

Next, we will formally analyze the condition that least-square loss fails to
identify the causal direction in the linear and nonlinear cases.

###  4.1 Linear Case

In the linear case, let $X\rightarrow Y$ denotes the causal direction from
cause $X$ to effect $Y$, and let $Y\rightarrow X$ denote the anti-causal
direction. Let $\sigma_{X}^{2}$ denote the variance of variable $X$. Let
$\beta_{Y|X}$ be the linear regression coefficient obtained from regressing
$Y$ on $X$, and let $\sigma_{Y|X}^{2}=\text{Var}(Y-\beta_{Y|X}\cdot X)$ be the
variance of the residual of regressing $Y$ on $X$.

Then, for the causal direction, we have:

$\begin{array}[]{c}\sigma_{X|\phi}^{2}=\sigma_{N_{X}}^{2},\quad\beta_{Y|X}=\frac{\alpha\sigma_{X}^{2}}{\sigma_{X}^{2}}=\alpha,\quad\sigma_{Y|X}^{2}=\sigma_{N_{Y}}^{2}.\end{array}$
(4)

For the reverse direction, we have:

$\displaystyle\sigma^{2}_{Y|\phi}$
$\displaystyle=\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}},\quad\beta_{X|Y}=\frac{\alpha\sigma^{2}_{N_{X}}}{\sigma^{2}_{Y|\phi}}=\frac{\alpha\sigma^{2}_{N_{X}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}},$
(5) $\displaystyle\sigma^{2}_{X|Y}$
$\displaystyle=\text{Var}(X-\beta_{X|Y}Y)=\frac{\sigma^{2}_{N_{X}}\sigma^{2}_{N_{Y}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}.$

Note that NOTEARS will prefer the direction that has a smaller least-square
loss. Thus, we are interested to see in what condition NOTEARS will fail,
i.e., the least-square loss satisfying $LS_{X\to Y}>LS_{Y\to X}$. By using the
above derivation, we can show such a condition in the following theorem.

######  Theorem 1.

Let $X\to Y$ be the causal direction following the data generation mechanism
$X=N_{X},Y=\alpha X+N_{Y}$. The least-square loss will fail to identify the
correct causal direction if the causal strength $\alpha$, the noise variances
$\sigma_{N_{X}}^{2}$ and $\sigma_{N_{Y}}^{2}$ satisfy the following
inequality:

$\alpha^{2}<1-\frac{\sigma_{N_{Y}}^{2}}{\sigma_{N_{X}}^{2}}.$ (6)

###### Proof.

Note that the least-square loss in the additive noise model is equivalent to
the summation of estimated noise variance. Thus, using the Eq. 4, the least-
square loss in the direction of $\displaystyle X\rightarrow Y$ is given as
follows:

$LS_{X\rightarrow
Y}=\mathbb{E}\left[(X-0)^{2}\right]+\mathbb{E}\left[(Y-\beta_{Y|X}\cdot
X)^{2}\right]=\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}},$ (7)

where $\displaystyle\beta_{Y|X}$ is the coefficient obtained from linear
regression. Similarly, using Eq. 5, we obtain

$LS_{Y\rightarrow X}=\mathbb{E}\left[(X-\beta_{X|Y}\cdot
Y)^{2}\right]+\mathbb{E}\left[(Y-0)^{2}\right]=\frac{\sigma^{2}_{N_{Y}}\sigma^{2}_{N_{X}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}+\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}.$
(8)

Using Eq. 7 and 8, we consider the condition that the loss in causal direction
is larger than the loss in the reverse direction, i.e., $LS_{X\rightarrow
Y}>LS_{Y\rightarrow X}$, and we obtain 6. The details of derivation are given
in Appendix A. ∎

Theorem 1 indicates that, under the linear additive noise model assumption,
the least-square loss could fail to identify the correct causal direction
unless the inequality 6 does not hold. It is also interesting to see that if
all noises follow the standard Gaussian distribution, we have
$\sigma_{N_{X}}^{2}=\sigma_{N_{Y}}^{2}$, then the causal direction can still
be correctly identified because the inequality $\alpha^{2}<0$ must not hold.
However, it is unrealistic that all noises follow the same standard Gaussian
distribution in the real world. And a simple standardization that is usually
used for preprocessing will invalidate the assumption. In addition, Theorem 1
also explains the experiments in Figure 1 that why the change of causal
strength or the noise variance will change the causal orientation.

Therefore, we conclude that the least-square loss is not a suitable score for
learning linear causal structure.

###  4.2 Nonlinear Case

In the nonlinear case, we study SCM between $X\to Y$ with the nonlinear
additive form $Y=f(X)+N_{Y}$. To show the limitation of least-square loss in a
nonlinear case, we aim to show that using the least-square loss can be viewed
as minimizing the mutual information between the noise and cause variables but
with bias such that the independence measure will be not reliable, leading to
identify the wrong direction. To do so, we first show that maximizing
likelihood is equivalent to minimizing the mutual information between the
cause and the noise variables (Lemma 1). Then, we further show that maximizing
the likelihood under the standard Gaussian noise assumption is equivalent to
minimizing the least-square loss (Lemma 2). Combing with Lemma 1 and Lemma 2,
we conclude that using least-square loss, the mutual information will be
biased (Lemma 3). Finally, by considering the bias in mutual information, we
provide the condition that least-square loss fails in Theorem 2.

To be clear, let $\hat{f}$ denote the function of regressing $Y$ on $X$, and
$\hat{g}$ denote the function of regressing $X$ on $Y$. Keep in mind that
$N_{X}$ and $N_{Y}$ are the noise variables. To avoid ambiguity, we denote
that: $\hat{N}_{Y}:=Y-\hat{f}(X)$ and $\hat{N}_{X}:=X-\hat{g}(Y)$. Let
$J_{X\to Y}$ denote the Jacobian matrix of the transformation from
$(X,N_{Y})^{T}$ to $(X,Y)^{T}$ , i.e., $J_{X\to
Y}=\begin{pmatrix}\frac{\partial X}{\partial X}&\frac{\partial X}{\partial
N_{Y}}\\\ \frac{\partial Y}{\partial X}&\frac{\partial Y}{\partial
N_{Y}}\end{pmatrix}$, and let $J_{Y\to X}$ denote the Jacobian matrix of the
transformation from $(Y,\hat{N}_{X})^{T}$ to $(Y,X)^{T}$.

######  Lemma 1 (Theorem 2 in Zhang et al. (2015)).

Given samples $\left\\{x^{(i)},y^{(i)}\right\\}^{m}_{i=1}$, and the causal
model $Y=f(X,N_{Y};\theta)$ with any parameters $\displaystyle\theta$ and the
estimated residual $\hat{N}_{Y}$, the log-likelihood $l_{X\rightarrow
Y}(\theta)$ and the mutual information between $X$ and $\hat{N}_{Y}$ are
related in the following way:

$I(X,\hat{N}_{Y};\theta)=\frac{1}{m}\sum^{m}_{i=1}\log
p(X=x^{(i)},Y=y^{(i)})-\frac{1}{m}l_{X\rightarrow Y}(\theta),$ (9)

where

$l_{X\rightarrow Y}(\theta)=\sum^{m}_{i=1}\log p(X=x^{(i)})+\sum^{m}_{i=1}\log
p(N_{Y}=\hat{n}^{(i)}_{Y};\theta)-\sum^{m}_{i=1}\log\left|\frac{\partial
f}{\partial N_{Y}}|_{N_{Y}=\hat{n}^{(i)}_{Y}}\right|.$ (10)

Lemma 1 builds a connection between the log-likelihood and the mutual
information. In the following Lemma, we further build the connection between
log-likelihood and the least-square loss.

######  Lemma 2.

For the additive noise model $Y=f(X)+N_{Y}$, maximizing the log-likelihood
$l_{X\to Y}(\theta)$ with the standard Gaussian noise assumption:

$\displaystyle l_{X\rightarrow
Y}(\theta)=\sum^{m}_{i=1}\log\left(\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(x^{(i)}-0)^{2}}{2}\right)\right)+\sum^{m}_{i=1}\log\left(\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(y^{(i)}-\hat{f}(x^{(i)};\theta))^{2}}{2}\right)\right)$
$\displaystyle-\sum^{m}_{i=1}\log\left|\frac{\partial f}{\partial
N_{Y}}|_{N_{Y}=\hat{n}^{(i)}_{Y}}\right|.$

is equivalent to minimizing the least-square loss $LS_{X\to Y}$:

$LS_{X\rightarrow
Y}=\mathbb{E}\left[(X-0)^{2}\right]+\mathbb{E}\left[(Y-\hat{f}(X))^{2}\right].$

Based on Lemma 1 and Lemma 2, one can straightforward develop the following
Lemma that using the least-square loss will introduce the standard Gaussian
distribution in the mutual information.

######  Lemma 3.

For the additive noise model, minimizing the least-square loss is equivalent
to minimizing the mutual information under the standard Gaussian noise
assumption with the following form:

$\displaystyle I_{q}(X,\hat{N}_{Y};\theta)$
$\displaystyle=\frac{1}{m}\left(\sum^{m}_{i=1}\log
p(X=x^{(i)},Y=y^{(i)})-\sum^{m}_{i=1}\log q(X=x^{(i)})-\sum^{m}_{i=1}\log
q(N_{Y}=\hat{n}^{(i)}_{Y};\theta)\right),$ (11)
$\displaystyle=\mathbb{E}_{x,y\sim p(x,y)}[\log p(X=x,Y=y)-\log q(X=x)-\log
q(N_{Y}=\hat{n}_{Y};\theta)].$

where $\displaystyle q$ is the density function of standard Gaussian
distribution.

Based on Lemma 3, in the following theorem, we show that in what condition the
least-square loss will fail.

######  Theorem 2.

Let $X\to Y$ be the causal direction following the data generation mechanism
$Y=f(X)+N_{Y}$. The causal direction is non-identifiable using least-square
loss if the following inequality holds:

$\begin{array}[]{l}-\int p(X)\log q(X)-\int p(N_{Y})\log q(N_{Y})>-\int
p(Y)\log q(Y)-\int p(\hat{N}_{X})\log q(\hat{N}_{X}),\end{array}$ (12)

where q is the density function of standard Gaussian distribution.

###### Proof.

Without loss of generality, we assume the identifiability condition of the
additive noise model holds. Then, we have $I(X,N_{Y})<I(Y,\hat{N}_{X})$.
However, based on Lemma 3, using the least-square loss, the mutual information
becomes:

$\displaystyle I_{q}(X,N_{Y})$ $\displaystyle=-\int p(X)\log q(X)-\int
p(N_{Y})\log q(N_{Y})+\int p(X,Y)\log p(X,Y),$ (13) $\displaystyle
I_{q}(Y,\hat{N}_{X})$ $\displaystyle=-\int p(Y)\log q(Y)-\int
p(\hat{N}_{X})\log q(\hat{N}_{X})+\int p(X,Y)\log p(X,Y),$

where $q\sim N(0,I)$. In this cases, the inequality $\displaystyle
I_{q}(X,N_{Y})<I_{q}(Y,\hat{N}_{X})$ is not necessary holds. In fact, by
solving the inequality $\displaystyle I_{q}(X,N_{Y})>I_{q}(Y,\hat{N}_{X})$
using 13, we obtain 12, which is the condition that least-square fails. ∎

Theorem 2 shows that if inequality 12 holds, the least-square loss could fail
especially when Kullback-Leibler divergence between distribution $p(X)$ or
$p(N_{Y})$ and the standard Gaussian distribution is large.

For example, let $Y=\text{sigmoid}(X)+N_{Y}$ where
$X\sim\text{Uniform}\left(-C,C\right)$ and
$N_{Y}\sim\text{Uniform}\left(-1,1\right)$. In this case, we can simply
increase $C$ such that $-\int p(X)\log q(X)$ will tend to infinity while $Y$
and $\hat{N}_{X}$ are bounded due to the sigmoid function, such that the right
hand side in 12 are also bound, and hence the inequality must hold.

Therefore, it is necessary to use the correct distribution setting $q=p$ to
obtain the correct result, which inspires a way that uses the entropy-based
loss instead.

##  5 Structure Learning Using Entropy-based Loss

As discussed in Section 4, using least-square loss, which is equivalent to
assuming the distribution of all noise terms is standard Gaussian
distribution, will result in the incorrect causal identification. The reason
for these errors is the mismatch of the truth distribution and the assumed
standard Gaussian distribution.

Instead of using the square operator in least-square loss, we can replace it
with the entropy such that:

$\min\sum_{i=1}^{d}H(N_{i})+\lambda\left\|W\right\|_{1}\quad\text{subject
to}\quad tr\left(e^{W\circ W}\right)-d=0,$ (14)

where $H\left(N_{i}\right)=-\int p\left(N_{i}\right)\log p\left(N_{i}\right)$
denotes the entropy of variable $N_{i}$.

To explain why we use entropy-based loss, we first show that the entropy-based
loss is consistent with the log-likelihood score under the additive noise
model (see Theorem 3). Based on Theorem 3, by utilizing identifiability
results in terms of likelihood, it is easy to show that the entropy-based loss
is also identifiable.

######  Theorem 3.

In additive noise model, the entropy-based score has a consistency with the
log-likelihood score when the sample sizes $m\rightarrow\infty$, that is

$\frac{1}{m}\sum^{m}_{j=1}\sum^{d}_{i=1}\log
p\left(x^{(j)}_{i}|x^{(j)}_{pa_{\mathcal{G}}(i)}\right)=-\sum^{d}_{i=1}H(N_{i}).$
(15)

Theorem 3 indicates that, under the additive noise model, the entropy-based
loss is consistent with the log-likelihood score, while Peters et al. (2014)
had proven that under the ANM model, the log-likelihood is able to distinguish
the causal direction under a mild assumption. Hence, the entropy-based loss
has the same property, which is illustrated in the following corollary:

######  Corollary 1.

For each pair of additive noise model $Y=f(X)+N_{Y}$, if
$v^{\prime\prime}(Y-f(X))f^{\prime}(X)\neq 0$ and the following condition
holds:

$\small\xi^{\prime\prime\prime}=\xi^{\prime\prime}\left(-\frac{v^{\prime\prime\prime}f^{\prime}}{v^{\prime\prime}}+\frac{f^{\prime\prime}}{f^{\prime}}\right)-2v^{\prime\prime}f^{\prime\prime}f^{\prime}+v^{\prime}f^{\prime\prime\prime}+\frac{v^{\prime}v^{\prime\prime\prime}f^{\prime\prime}f^{\prime}}{v^{\prime\prime}}-\frac{v^{\prime}(f^{\prime\prime})^{2}}{f^{\prime}},$
(16)

where $v\coloneq\log p(N_{Y})$, $\xi\coloneq\log p(X)$, then using the
entropy-based loss each pair of additive model is identifiable and the
following inequality holds:

$H\left(X\right)+H\left(\hat{N}_{Y}\right)<H\left(Y\right)+H\left(\hat{N}_{X}\right).$
(17)

Corollary 1 illustrates the validity of the entropy-based loss under the
additive noise model, which is based on the consistency with the likelihood
score, indicating that by comparing the entropy of residuals between causal
and anti-causal direction, we can identify the cause and effect variables in
most cases.

##  6 Experiments

To investigate the effectiveness of the entropy loss, we compare it with
baseline algorithms on both synthetic data and real-world data. The synthetic
data are generated from linear non-Gaussian data and nonlinear additive noise
data, respectively. The following algorithms that use the least-square loss
for data reconstruction are taken as the baseline: NOTEARS (Zheng et al.,
2018), DAG-GNN (Yu et al., 2019), NOTEARS-MLP (Zheng et al., 2020). We reuse
the parameter settings for those baseline algorithms in their original papers
and codes.

For the differential entropy estimator, it turn out that the method proposed
by Hyvärinen (1998) is well suit for this problem. in which the entropy can be
estimated using following method:

$H(X)\approx
H(\nu)-\left[k_{1}\left(E\left\\{\bar{G}_{1}(X)\right\\}\right)^{2}+k_{2}\left(E\left\\{\bar{G}_{2}(X)\right\\}-E\left\\{\bar{G}_{2}(\nu)\right\\}\right)^{2}\right],$
(18)

where $k_{1}=36/\left(8\sqrt{3}-9\right)$,
$k_{2}=24/\left(16\sqrt{3}-27\right)$,
$\bar{G}_{1}(X)=X\exp\left(-X^{2}/2\right)$,
$\bar{G}_{2}(X)=\exp\left(-X^{2}/2\right)$, $\bar{G}_{2}(\nu)=\sqrt{1/2}$, and
$H(\nu)=\frac{1}{2}\left(1+\log(2\pi)\right)$.

Our implementation, denoted by ours (linear) and ours-MLP (nonlinear), is
based on the codes of NOTEARS (Zheng et al., 2018) and NOTEARS-MLP (Zheng et
al., 2020), respectively. Following the previous works, the numerical optimize
algorithms L-BFGS (Byrd et al., 1995) is used for optimization, and we will
prune the edges after training with a small threshold $\omega=0.3$ to rule out
cycle-inducing edges.

In linear non-Gaussian case, the data are generated according to the following
linear structure equations, $X_{i}=\Sigma_{j\in
pa_{\mathcal{G}}(i)}\beta_{ij}X_{j}+N_{j}$ with random causal strength
$\beta_{ij}\sim\text{Uniform}(-0.4,-0.8)\cup\text{Uniform}(0.4,0.8)$ and the
noise $N_{j}\sim\sigma_{N_{j}}*\text{Uniform}(-\sqrt{3},\sqrt{3})$ or
$N_{j}\sim\sigma_{N_{j}}*\text{Gumbel}(0,\frac{\sqrt{6}}{\pi})$ with random
standard deviation $\sigma_{N_{j}}\sim\text{Uniform}(0.5,1.0)$. In nonlinear
case, the data are generated according to the following nonlinear SCM:
$X_{i}=\text{tanh}\left(\Sigma_{j\in
pa_{\mathcal{G}}(i)}\beta_{1,ij}X_{j}\right)+\text{cos}\left(\Sigma_{j\in
pa_{\mathcal{G}}(i)}\beta_{2,ij}X_{j}\right)+\text{sin}\left(\Sigma_{j\in
pa_{\mathcal{G}}(i)}\beta_{3,ij}X_{j}\right)+N_{j}$ with random causal
strength $\beta_{k,ij}\sim\text{Uniform}(0.5,2.0)$ for $k=1,2,3$ and the noise
$N_{j}\sim\sigma_{N_{j}}*\text{Uniform}(-\sqrt{3},\sqrt{3})$ with random
standard deviation $\sigma_{N_{j}}\sim\text{Uniform}(0.5,1.0)$.

Following previous works, Structural Hamming Distance (SHD), False Discovery
Rate (FDR), and True Positive Rate (TPR) are recorded as the evaluation
metrics for all algorithms.

###  6.1 Synthetic Structures

In this section, we design a series of controlled experiments with respect to
the sample sizes and the number of variables on the synthetic random causal
structures. At each experiment, we will control one of the parameters while
fixing others. The range of the above controlled parameters are as follows:
the number of samples=$\\{200,400,\textbf{600},800,1000\\}$ and the number of
variables$=\\{5,10,\textbf{15},20,25\\}$ with in-degree=2 where the default
setting is marked as bold. All experiments will run at least ten times.

####  6.1.1 Linear Case

(a) Uniform noises with variables=15

(b) Uniform noises with samples=600

(c) Gumbel noises with samples=600

Figure 2: Structure recovery in linear data in terms of SHD, FDR and TPR to
the true graph: in Figure 2(a), the noises are Uniform, the number of
variables and edges are 15 and 30; in Figure 2(b), the noises are Uniform, and
the sample sizes are 600; in Figure 2(c), the noises are Gumbel, and the
sample sizes are 600.

In this part, we perform linear experiments on the synthetic linear non-
Gaussian data. Here, we test our method on Uniform and the Gumbel
distribution, respectively. Each algorithm runs 30 times for all experiments.

As shown in Figure 2, our method outperforms all baseline methods on the three
controlled experiments. Moreover, by using the entropy-based method, we have
the least deviation error. It is also interesting to see the significant
difference between ours and the NOTEARS, which verifies that the algorithm
could learn the reverse causal direction using least-square loss. In detail,
in Figure 2(a), as the sample size increase, all methods has the better
performance, and we can see that when the sample size is small, we still
outperform the other methods, which shows the robustness of our method. In
Figure 2(b) and 2(b), in general, the performance decrease as the number of
variables increases, but compared with the baseline methods, our method
decreases slowly. In addition, we can see that under different distributions,
our method has similar performance within experimental errors, which verifies
the Theorem 1.

####  6.1.2 Nonlinear Case

(a) Uniform noises with variables=15

(b) Uniform noises with samples=600

(c) Uniform noises with variables=15

Figure 3: Structure recovery in nonlinear data in terms of SHD, FDR and TPR to
the true graph: in figure 3(a), the noises are uniform with variance=3.0 and
the number of variables and edges are 15 and 30; in Figure 3(b), the noises
are Uniform with variance=3.0 and the number of sample sizes is 600; in Figure
3(c), the noises are Uniform, the number of sample sizes is 600, and the
number of variables and edges are 15 and 30.

In this part, we test our entropy-based loss in nonlinear data. To further
verify the Theorem 3 we generate the nonlinear data using Uniform distribution
with different variances. That is, the higher the variance, the more non-
Gaussianity. In the controlled experiments we further test our method on
different noise variance=$\\{1,2,\textbf{3},4,5\\}$. Each algorithm will run
at least ten times.

The results are given in Figure 3. Overall, our method generally outperforms
the baseline methods showing the robustness of our method in variable-varying
and variance-varying cases. In detail, in Figure 3(a), we notice that
comparing with the linear data, all methods require more data to acquire a
decent performance, but our method still outperforms the baseline methods. For
our method, 600 is a decent choice for the sample sizes. In figure 3(b), with
the growing number of variables, the performance of NOTEARS-MLP decreases
rapidly while our method remains stable. The reason is that as the number of
variables grows, the edges increase simultaneously, and the incorrectly
identified edges will also increase. In Figure 3(b), we can see that our
method is not sensitive to the noise variances while the performance of other
methods decreases rapidly as the variances grow. The reason is that the
variance controls the noise non-Gaussianity making it far from the standard
Gaussian. Therefore, the entropy-based loss is stable and identifiable under
any distribution while other methods will lose its identifiability.

###  6.2 Real-World Data

Following previous works, we use the real-world data set provided by Sachs et
al. (2005). This data set contains 7466 groups of continuous expression levels
of 20 measured molecules in human immune system cells. The consensus network,
accepted by the biological community, is represented by a Bayesian network
with 11 variables and 20 edges. Compared with NOTEARS that predicts 16 edges
with an SHD of 22, ours predicts 19 edges with an SHD of 21, which shows the
effectiveness of our method with higher recall and precision. Compared with
NOTEARS-MLP predicting 13 edges with an SHD of 16, and DAG-GNN predicting 18
edges with an SHD of 19, ours-MLP has a similar result predicting 13 edges
with an SHD of 16.

##  7 Conclusion

In this work, we re-examine the least-square loss for learning causal
structures in causal discovery and advocate to use a more general entropy-
based loss as the objective function. Our analysis shows that the least-square
loss is problematic to learn the underlying causal structure if the standard
Gaussian noise assumption does not hold. On the contrary, the entropy-based
loss is consistent with the log-likelihood score, which is robust and
effective for learning causal structures. Our experimental results verify the
proposed theoretical results and show the effectiveness of the proposed
entropy-based loss. A clear next step is to generalize the theory and
algorithms to a more general causal mechanism.

## References

  * Andersson et al. [1997] Steen A Andersson, David Madigan, Michael D Perlman, et al.  A characterization of markov equivalence classes for acyclic digraphs.  _Annals of statistics_ , 25(2):505–541, 1997\.
  * Blöbaum et al. [2018] Patrick Blöbaum, Dominik Janzing, Takashi Washio, Shohei Shimizu, and Bernhard Schölkopf.  Cause-effect inference by comparing regression errors.  In _International Conference on Artificial Intelligence and Statistics_ , pages 900–909. PMLR, 2018.
  * Blöbaum et al. [2019] Patrick Blöbaum, Dominik Janzing, Takashi Washio, Shohei Shimizu, and Bernhard Schölkopf.  Analysis of cause-effect inference by comparing regression errors.  _PeerJ Computer Science_ , 5:e169, 2019.
  * Byrd et al. [1995] Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A limited memory algorithm for bound constrained optimization.  _SIAM Journal on scientific computing_ , 16(5):1190–1208, 1995.
  * Cai et al. [2016] Ruichu Cai, Zhenjie Zhang, Zhifeng Hao, and Marianne Winslett.  Understanding social causalities behind human action sequences.  _IEEE transactions on neural networks and learning systems_ , 28(8):1801–1813, 2016.
  * Didelez and Pigeot [2001] Vanessa Didelez and Iris Pigeot.  Judea pearl: Causality: Models, reasoning, and inference.  _Politische Vierteljahresschrift_ , 42(2):313–315, 2001.
  * Ghysels et al. [2016] Eric Ghysels, Jonathan B Hill, and Kaiji Motegi.  Testing for granger causality with mixed frequency data.  _Journal of Econometrics_ , 192(1):207–230, 2016\.
  * Grosse-Wentrup et al. [2016] Moritz Grosse-Wentrup, Dominik Janzing, Markus Siegel, and Bernhard Schölkopf.  Identification of causal relations in neuroimaging data with latent confounders: An instrumental variable approach.  _NeuroImage_ , 125:825–833, 2016.
  * Hoyer et al. [2008] Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf.  Nonlinear causal discovery with additive noise models.  _Advances in neural information processing systems_ , 21:689–696, 2008.
  * Hyvärinen [1998] Aapo Hyvärinen.  New approximations of differential entropy for independent component analysis and projection pursuit.  In _Proceedings of the 1997 Conference on Advances in Neural Information Processing Systems 10_ , NIPS ’97, page 273–279, Cambridge, MA, USA, 1998. MIT Press.  ISBN 0262100762.
  * Lachapelle et al. [2019] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien.  Gradient-based neural dag learning.  _arXiv preprint arXiv:1906.02226_ , 2019.
  * Ng et al. [2019a] Ignavier Ng, Zhuangyan Fang, Shengyu Zhu, Zhitang Chen, and Jun Wang.  Masked gradient-based causal structure learning.  _arXiv preprint arXiv:1910.08527_ , 2019a.
  * Ng et al. [2019b] Ignavier Ng, Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang.  A graph autoencoder approach to causal structure learning.  _arXiv preprint arXiv:1911.07420_ , 2019b.
  * Peters and Bühlmann [2014] Jonas Peters and Peter Bühlmann.  Identifiability of gaussian structural equation models with equal error variances.  _Biometrika_ , 101(1):219–228, 2014.
  * Peters et al. [2014] Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf.  Causal discovery with continuous additive noise models.  _Journal of Machine Learning Research_ , 15(58):2009–2053, 2014.  URL http://jmlr.org/papers/v15/peters14a.html.
  * Sachs et al. [2005] Karen Sachs, Omar Perez, Dana Pe’er, Douglas A Lauffenburger, and Garry P Nolan.  Causal protein-signaling networks derived from multiparameter single-cell data.  _Science_ , 308(5721):523–529, 2005.
  * Shimizu et al. [2006] Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan.  A linear non-gaussian acyclic model for causal discovery.  _Journal of Machine Learning Research_ , 7(10), 2006.
  * Spirtes et al. [2000] Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman.  _Causation, prediction, and search_.  MIT press, 2000.
  * Yu et al. [2019] Yue Yu, Jie Chen, Tian Gao, and Mo Yu.  Dag-gnn: Dag structure learning with graph neural networks.  In _International Conference on Machine Learning_ , pages 7154–7163. PMLR, 2019.
  * Zhang and Hyvärinen [2009] Kun Zhang and Aapo Hyvärinen.  On the identifiability of the post-nonlinear causal model.  In Jeff A. Bilmes and Andrew Y. Ng, editors, _UAI 2009, Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence, Montreal, QC, Canada, June 18-21, 2009_ , pages 647–655. AUAI Press, 2009.
  * Zhang et al. [2015] Kun Zhang, Zhikun Wang, Jiji Zhang, and Bernhard Schölkopf.  On estimation of functional causal models: general results and application to the post-nonlinear causal model.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ , 7(2):1–22, 2015.
  * Zheng et al. [2018] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing.  Dags with no tears: Continuous optimization for structure learning.  _Advances in Neural Information Processing Systems_ , 31:9472–9483, 2018.
  * Zheng et al. [2020] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing.  Learning sparse nonparametric dags.  In _International Conference on Artificial Intelligence and Statistics_ , pages 3414–3425. PMLR, 2020.

## Supplementary Material

##  Appendix A Proof of Theorem 1

######  Theorem 1.

Let $X\to Y$ be the causal direction following the data generation mechanism
$X=N_{X},Y=\alpha X+N_{Y}$. The least-square loss will fail to identify the
correct causal direction if the causal strength $\alpha$, the noise variances
$\sigma_{N_{X}}^{2}$ and $\sigma_{N_{Y}}^{2}$ satisfy the following
inequality:

$\alpha^{2}<1-\frac{\sigma_{N_{Y}}^{2}}{\sigma_{N_{X}}^{2}}.$ (A.1)

###### Proof.

Then, for the causal direction, we have:

$\begin{array}[]{c}\sigma_{X|\phi}^{2}=\sigma_{N_{X}}^{2},\quad\beta_{Y|X}=\frac{\alpha\sigma_{X}^{2}}{\sigma_{X}^{2}}=\alpha,\quad\sigma_{Y|X}^{2}=\sigma_{N_{Y}}^{2}.\end{array}$
(A.2)

For the anti-causal direction, we have

$\displaystyle\sigma^{2}_{Y|\phi}$
$\displaystyle=\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}},\quad\beta_{X|Y}=\frac{\alpha\sigma^{2}_{N_{X}}}{\sigma^{2}_{Y|\phi}}=\frac{\alpha\sigma^{2}_{N_{X}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}$
(A.3) $\displaystyle\sigma^{2}_{X|Y}$
$\displaystyle=\text{Var}(X-\beta_{X|Y}Y)=\frac{\sigma^{2}_{N_{X}}\sigma^{2}_{N_{Y}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}.$

Note that the least-square loss in the additive noise model is equivalent to
the summation of estimated noise variance. Thus, using the Eq. A.2, the least-
square loss in the direction of $\displaystyle X\rightarrow Y$ is given as
follows:

$LS_{X\rightarrow
Y}=\mathbb{E}\left[(X-0)^{2}\right]+\mathbb{E}\left[(Y-\beta_{Y|X}\cdot
X)^{2}\right]=\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}},$ (A.4)

where $\displaystyle\beta_{Y|X}$ is the coefficient obtained from linear
regression. Similarly, using Eq. A.3, we obtain

$LS_{Y\rightarrow X}=\mathbb{E}\left[(X-\beta_{X|Y}\cdot
Y)^{2}\right]+\mathbb{E}\left[(Y-0)^{2}\right]=\frac{\sigma^{2}_{N_{Y}}\sigma^{2}_{N_{X}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}+\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}.$
(A.5)

Using Eq. A.4 and A.5, we consider the condition that the loss in causal
direction is larger than the loss in the reverse direction, i.e.,
$LS_{X\rightarrow Y}>LS_{Y\rightarrow X}$, and we have the following
inequality:

$\displaystyle\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}$
$\displaystyle>\frac{\sigma^{2}_{N_{Y}}\sigma^{2}_{N_{X}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}+\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}$
(A.6) $\displaystyle\sigma^{2}_{N_{X}}$
$\displaystyle>\left(\frac{\sigma^{2}_{N_{Y}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}+\alpha^{2}\right)\sigma^{2}_{N_{X}}$
$\displaystyle\alpha^{2}$
$\displaystyle<1-\frac{\sigma^{2}_{N_{Y}}}{\alpha^{2}\sigma^{2}_{N_{X}}+\sigma^{2}_{N_{Y}}}$
$\displaystyle\alpha^{2}$
$\displaystyle<\frac{\alpha^{2}\sigma_{N_{X}}^{2}+\sigma_{N_{Y}}^{2}-\sigma_{N_{Y}}^{2}}{\alpha^{2}\sigma_{N_{X}}^{2}+\sigma_{N_{Y}}^{2}}$
$\displaystyle\alpha^{2}$
$\displaystyle<\frac{\alpha^{2}}{\alpha^{2}+\frac{\sigma_{N_{Y}}^{2}}{\sigma_{N_{X}}^{2}}}$
$\displaystyle\alpha^{2}$
$\displaystyle<1-\frac{\sigma_{N_{Y}}^{2}}{\sigma_{N_{X}}^{2}}.$

∎

##  Appendix B Proof of Lemma 2

######  Lemma 1 (Theorem 2 in Zhang et al. [2015]).

Given samples $\left\\{x^{(i)},y^{(i)}\right\\}^{m}_{i=1}$, and the causal
model $Y=f(X,N_{Y};\theta)$ with any parameters $\displaystyle\theta$ and the
estimated residual $\hat{N}_{Y}$, the log-likelihood $l_{X\rightarrow
Y}(\theta)$ and the mutual information between $X$ and $\hat{N}_{Y}$ are
related in the following way:

$I(X,\hat{N}_{Y};\theta)=\frac{1}{m}\sum^{m}_{i=1}\log
p(X=x^{(i)},Y=y^{(i)})-\frac{1}{m}l_{X\rightarrow Y}(\theta),$ (B.1)

where

$l_{X\rightarrow Y}(\theta)=\sum^{m}_{i=1}\log p(X=x^{(i)})+\sum^{m}_{i=1}\log
p(N_{Y}=\hat{n}^{(i)}_{Y};\theta)-\sum^{m}_{i=1}\log\left|\frac{\partial
f}{\partial N_{Y}}|_{N_{Y}=\hat{n}^{(i)}_{Y}}\right|.$ (B.2)

Lemma 1 builds a connection between the log-likelihood and the mutual
information. In the following Lemma, we further build the connection between
log-likelihood and the least-square loss.

######  Lemma 2.

For the additive noise model $Y=f(X)+N_{Y}$, maximizing the log-likelihood
$l_{X\to Y}(\theta)$ with the standard Gaussian noise assumption is equivalent
to minimizing the least-square loss.

###### Proof.

Based on Lemma 1, the log-likelihood under the standard Gaussian distribution
has the following form:

$\displaystyle l_{X\rightarrow
Y}(\theta)=\sum^{m}_{i=1}\log\left(\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(x^{(i)}-0)^{2}}{2}\right)\right)+\sum^{m}_{i=1}\log\left(\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(y^{(i)}-\hat{f}(x^{(i)};\theta))^{2}}{2}\right)\right)$
(B.3) $\displaystyle-\sum^{m}_{i=1}\log\left|\frac{\partial f}{\partial
N_{Y}}|_{N_{Y}=\hat{n}^{(i)}_{Y}}\right|.$

Note that, under the additive noise model, we have
$\displaystyle\frac{\partial f}{\partial N_{Y}}=1$, and thus last term in B.3
will be vanish. Then, log-likelihood can be rewritten as follows:

$l_{X\rightarrow
Y}(\theta)=-\frac{1}{2}\sum^{m}_{i=1}\left[(x^{(i)}-0)^{2}+((y^{(i)}-\hat{f}(x^{(i)};\theta))^{2}\right]+const.$
(B.4)

Comparing with the least-square loss $LS_{X\rightarrow Y}$ that used in
NOTEARS:

$LS_{X\rightarrow
Y}=\mathbb{E}\left[(X-0)^{2}\right]+\mathbb{E}\left[(Y-\hat{f}(X))^{2}\right],$
(B.5)

we can see that maximizing the log-likelihood with standard Gaussian noise
assumption is equivalent to minimizing the least-square loss. ∎

##  Appendix C Proof of Lemma 3

######  Lemma 3.

For the additive noise model, minimizing the least-square loss is equivalent
to minimizing the mutual information under the standard Gaussian noise
assumption with the following form:

$\displaystyle I_{q}(X,\hat{N}_{Y};\theta)$
$\displaystyle=\frac{1}{m}\left(\sum^{m}_{i=1}\log
p(X=x^{(i)},Y=y^{(i)})-\sum^{m}_{i=1}\log q(X=x^{(i)})-\sum^{m}_{i=1}\log
q(N_{Y}=\hat{n}^{(i)}_{Y};\theta)\right),$ (C.1)
$\displaystyle=\mathbb{E}_{x,y\sim p(x,y)}[\log p(X=x,Y=y)-\log q(X=x)-\log
q(N_{Y}=\hat{n}_{Y};\theta)].$

where $\displaystyle q$ is the density function of standard Gaussian
distribution.

###### Proof.

Based on Lemma 2, as Eq. B.3 shown, using the least-square loss, the
distribution in the log-likelihood will become Gaussian. Then, by substituting
B.3 into B.1 we obtain C.1 which is the mutual information under the standard
Gaussian noise assumption. ∎

##  Appendix D Proof of Theorem 3

######  Theorem 3.

In additive noise model, the entropy-based score has a consistency with the
log-likelihood score when the sample sizes $m\rightarrow\infty$, that is

$\frac{1}{m}\sum^{m}_{j=1}\sum^{d}_{i=1}\log
p\left(x^{(j)}_{i}|x^{(j)}_{pa_{\mathcal{G}}(i)}\right)=-\sum^{d}_{i=1}H(N_{i}).$
(D.1)

###### Proof.

Based on Lemma 1 we generalize the log-likelihood from two variables case to
the multivariate case as follows:

$\displaystyle\frac{1}{m}\sum^{m}_{j=1}\sum^{d}_{i=1}\log
p\left(x^{(j)}_{i}|x^{(j)}_{pa_{\mathcal{G}}(i)}\right)$
$\displaystyle=\frac{1}{m}\sum^{m}_{j=1}\sum^{d}_{i=1}\log
p(N_{i}=\hat{n}^{(j)}_{i})-\frac{1}{m}\sum^{m}_{j=1}\sum^{d}_{i=1}\log\left|\frac{\partial
f_{i}}{\partial N_{i}}|_{N_{i}=\hat{n}^{(j)}_{i}}\right|$ (D.2)
$\displaystyle=\frac{1}{m}\sum^{m}_{j=1}\sum^{d}_{i=1}\log
p\left(N_{i}=\hat{n}^{(j)}_{i}\right)=\sum^{d}_{i=1}-H(N_{i}).$

The second equality is based on the property of the additive noise model in
which $\frac{\partial f_{i}}{\partial N_{i}}=1$. Then, as
$m\rightarrow\infty$, we obtain the last equality. ∎

##  Appendix E Proof of Corollary 1

######  Lemma 4 (Theorem 1 in Hoyer et al. [2008]).

For the additive noise model $Y=f(X)+N_{Y}$, there is a the forward model of
the form $p(X,Y)=p(Y-f(X))p(X)$. If there is a backward model of the same form
$p(X,Y)=p(X-f(Y))p(Y)$, then for all $X$, $Y$ with the three time
differentiable $f$ and $v^{\prime\prime}(Y-f(X))f^{\prime}(X)\neq 0$, the
following equality holds:

$\xi^{\prime\prime\prime}=\xi^{\prime\prime}\left(-\frac{v^{\prime\prime\prime}f^{\prime}}{v^{\prime\prime}}+\frac{f^{\prime\prime}}{f^{\prime}}\right)-2v^{\prime\prime}f^{\prime\prime}f^{\prime}+v^{\prime}f^{\prime\prime\prime}+\frac{v^{\prime}v^{\prime\prime\prime}f^{\prime\prime}f^{\prime}}{v^{\prime\prime}}-\frac{v^{\prime}(f^{\prime\prime})^{2}}{f^{\prime}},$
(E.1)

where $v\coloneq\log p(N_{Y})$, $\xi\coloneq\log p(X)$.

Lemma 4 guarantees that the asymmetry of independence only vanishes under
strictly conditions. In most cases, such asymmetry exists and the ANM can be
identified. Based on Lemma 4, we prove the asymmetry of entropy also exists
for most cases.

######  Corollary 1.

For each pair of additive noise model $Y=f(X)+N_{Y}$, if
$v^{\prime\prime}(Y-f(X))f^{\prime}(X)\neq 0$ and the following condition
holds:

$\xi^{\prime\prime\prime}=\xi^{\prime\prime}\left(-\frac{v^{\prime\prime\prime}f^{\prime}}{v^{\prime\prime}}+\frac{f^{\prime\prime}}{f^{\prime}}\right)-2v^{\prime\prime}f^{\prime\prime}f^{\prime}+v^{\prime}f^{\prime\prime\prime}+\frac{v^{\prime}v^{\prime\prime\prime}f^{\prime\prime}f^{\prime}}{v^{\prime\prime}}-\frac{v^{\prime}(f^{\prime\prime})^{2}}{f^{\prime}},$
(E.2)

where $v\coloneq\log p(N_{Y})$, $\xi\coloneq\log p(X)$, then using the
entropy-based loss each pair of additive model is identifiable and has the
form
$H\left(X\right)+H\left(\hat{N}_{Y}\right)<H\left(Y\right)+H\left(\hat{N}_{X}\right).$

###### Proof.

To be clear, let $\hat{f}$ denote the function of regressing Y on X with a
residual $\hat{N}_{Y}$, and let $\hat{g}$ denote the function of regressing X
on Y with a residual $\hat{N}_{X}$.

To prove corollary 1, we will compare the mutual information between
hypothetical cause and the regression residual in causal and anti-causal
direction. Assuming the backward model does not exist based on Lemma 4, the
noise and the hypothetical cause is independent, i.e., $N_{i}\Vbar
X_{pa_{\mathcal{G}}(i)}$ in the causal direction, but not vice versa. Such a
property can be represented by mutual information form as follows:

$\displaystyle I(X,\hat{N}_{Y})$ $\displaystyle<I(Y,\hat{N}_{X})$ (E.3)
$\displaystyle\int
p(X,\hat{N}_{Y})\log\frac{p(X,\hat{N}_{Y})}{p(X)p(\hat{N}_{Y})}$
$\displaystyle<\int
p(X,\hat{N}_{Y})\log\frac{p(Y,\hat{N}_{X})}{p(Y)p(\hat{N}_{X})}$
$\displaystyle-\mathbb{E}\log p(X)-\mathbb{E}\log
p(\hat{N}_{Y})+\mathbb{E}\log p(X,\hat{N}_{Y})$ $\displaystyle<-\mathbb{E}\log
p(Y)-\mathbb{E}\log p(\hat{N}_{X})+\mathbb{E}\log p(Y,\hat{N}_{X})$
$\displaystyle H\left(X\right)+H\left(\hat{N}_{Y}\right)+\mathbb{E}\log
p(X,\hat{N}_{Y})$
$\displaystyle<H\left(Y\right)+H\left(\hat{N}_{X}\right)+\mathbb{E}\log
p(Y,\hat{N}_{X})$

In fact, $p(X,\hat{N}_{Y})$ and $p(Y,\hat{N}_{X})$ can be transformed to
$p(X,Y)$, and then be eliminated in both sides. We also denote $X^{\prime}$
and $Y^{\prime}$ for conveniently and clearly expressing the transformations
as follows:

$\text{causal direction:}\left\\{\begin{aligned} X^{\prime}&=X,\\\
\hat{N}_{Y}&=Y-\hat{f}(X).\end{aligned}\right.\quad\quad\text{anti-causal
direction:}\left\\{\begin{aligned} Y^{\prime}&=Y,\\\
\hat{N}_{X}&=X-\hat{g}(Y).\end{aligned}\right.$ (E.4)

Then, the distribution transformations from $(X,Y)^{T}$ to
$(X^{\prime},\hat{N}_{Y})^{T}$ and $(X,Y)^{T}$ to
$(Y^{\prime},\hat{N}_{X})^{T}$ are as follows:

$p(X,Y)=p(X^{\prime},\hat{N}_{Y})\left|\det(J_{X\rightarrow
Y})\right|\quad\quad
p(X,Y)=p(Y^{\prime},\hat{N}_{X})\left|\det(J_{Y\rightarrow X})\right|$

and

$\displaystyle\ \left|\det(J_{X\rightarrow
Y})\right|=\left|\det\begin{pmatrix}\frac{\partial X^{\prime}}{\partial
X}&\frac{\partial X^{\prime}}{\partial Y}\\\
\frac{\partial\hat{N}_{Y}}{\partial X}&\frac{\partial\hat{N}_{Y}}{\partial
Y}\end{pmatrix}\right|=\left|\det\begin{pmatrix}\frac{\partial X}{\partial
X}&\frac{\partial X}{\partial Y}\\\ \frac{\partial(Y-\hat{f}(X))}{\partial
X}&\frac{\partial(Y-\hat{f}(X))}{\partial
Y}\end{pmatrix}\right|=\left|\det\begin{pmatrix}1&0\\\
-\frac{\partial\hat{f}(X)}{\partial X}&1\end{pmatrix}\right|=1,$ (E.5)
$\displaystyle\ \left|\det(J_{Y\rightarrow
X})\right|=\left|\det\begin{pmatrix}\frac{\partial Y^{\prime}}{\partial
X}&\frac{\partial Y^{\prime}}{\partial Y}\\\
\frac{\partial\hat{N}_{X}}{\partial X}&\frac{\partial\hat{N}_{X}}{\partial
Y}\end{pmatrix}\right|=\left|\det\begin{pmatrix}\frac{\partial Y}{\partial
X}&\frac{\partial Y}{\partial Y}\\\ \frac{\partial(X-\hat{g}(Y))}{\partial
X}&\frac{\partial(X-\hat{g}(Y))}{\partial
Y}\end{pmatrix}\right|=\left|\det\begin{pmatrix}0&1\\\
1&-\frac{\partial\hat{g}(Y)}{\partial Y}\end{pmatrix}\right|=1,$

where $J_{X\rightarrow Y}$ denotes the Jacobian matrix of the transformation
from $(X,Y)^{T}$ to $(X^{\prime},\hat{N}_{Y})^{T}$ and $J_{Y\rightarrow X}$
vice versa.

Thus using Eq. E.5, we rewrite the last terms in both sides of E.3:

$\displaystyle\mathbb{E}\log p(X,\hat{N}_{Y})=\mathbb{E}\log
p(X^{\prime},\hat{N}_{Y})=\mathbb{E}(\log p(X,Y)/\left|\det(J_{X\to
Y})\right|)=\mathbb{E}\log p(X,Y)$ (E.6) $\displaystyle\mathbb{E}\log
p(Y,\hat{N}_{X})=\mathbb{E}\log p(Y^{\prime},\hat{N}_{X})=\mathbb{E}(\log
p(X,Y)/\left|\det(J_{Y\to X})\right|)=\mathbb{E}\log p(X,Y)$

Using E.3 and E.6, we obtain:

$H\left(X\right)+H\left(\hat{N}_{Y}\right)<H\left(Y\right)+H\left(\hat{N}_{X}\right),$
(E.7)

which finishes the proof. ∎

##  Appendix F Additional Experimental Details

###  F.1 Algorithm Parameter Settings

Table 1: Parameter Settings Parameter | Symbol | Value | Applicable to
---|---|---|---
threshold on $W$ | $\omega$ | 0.3 | all
acyclity penalty | $h(W)$ | $tr\left(e^{W\circ W}\right)-d$ | NOTEARS, NOTEARS-MLP ,
ours, ours-MLP
acyclity penalty | $h(W)$ | $tr\left((I+\frac{W\circ W}{d})^{d}\right)-d$ | DAG-GNN
$h$ tolerance | $\epsilon$ | $10^{-8}$ | all
$h$ progress rate | $c$ | $0.25$ | all
initial Lagrange multiplier | $\alpha_{0}$ | $0$ | all
$\rho$ increase factor |  | 10 | all
$\rho$ maximum |  | $10^{16}$ | all
optimize algorithms |  | L-BFGS | NOTEARS, NOTEARS-MLP ,
ours, ours-MLP
optimize algorithms |  | Adam | DAG-GNN
seed |  | $123$ | all

Parameter settings for all algorithms are shown in Table 1. We use the least-
squares loss
$\frac{1}{2m}\|\mathbf{X}-\mathbf{X}W\|^{2}_{F}+\lambda\left\|W\right\|_{1}$
regardless of the noise type for NOTEARS [Zheng et al., 2018] and NOTEARS-MLP
[Zheng et al., 2020]. We use the negative ELBO as the object function under
standard Gaussian noise assumption for DAG-GNN [Yu et al., 2019]. For our and
our-MLP, we use entropy-based loss as the objective function. The NOTEARS code
is available at https://github.com/xunzheng/notears and the DAG-GNN code is
available at https://github.com/fishmoon1234/DAG-GNN . Our code is available
at Supplemental Material.

###  F.2 Computing Environment

Solution times were obtained using a single CPU (Intel(R) Xeon(R) CPU E5-2620
v3 @ 2.40GHz).

###  F.3 License

We strictly follow the agreement of Apache-2.0 license used by NOTEARS code
and do not use NOTEARS code commercially.


# On conformal points of area preserving maps and related topics

Peter Albers111 Mathematisches Institut, Universität Heidelberg, 69120
Heidelberg, Germany<EMAIL_ADDRESS>Serge Tabachnikov222
Department of Mathematics, Pennsylvania State University, University Park, PA
16802, USA<EMAIL_ADDRESS>
###### Abstract

In this article we consider area preserving diffeomorphisms of planar domains,
and we are interested in their conformal points, i.e., points at which the
derivative is a similarity. We present some conditions that guarantee
existence of conformal points for the infinitesimal problem of Hamiltonian
vector fields as well as for what we call moderate symplectomorphisms of
simply connected domains. We also link this problem to the Carathéodory and
Loewner conjectures.

##  1 Intoduction

A conformal point of a diffeomorphism $F:{\mathbb{R}}^{2}\to{\mathbb{R}}^{2}$
is a point at which the derivative $dF$ is a similarity. In this article we
consider area preserving diffeomorphisms of planar domains, and we are
interested in their conformal points. Although such an area preserving
diffeomorphism may be free of conformal points (see Example 2.1 below), we
present some conditions on $F$ that guarantee their existence.

There are several motivations for this study. It was sparked by a question
posed on MathOverflow [20], apparently motivated by the elasticity theory: Is
it true that an area preserving diffeomorphism of the closed unit disc
$D^{2}\subset{\mathbb{R}}^{2}$ possesses conformal points?

We start with the infinitesimal version of the problem where a
symplectomorphism is replaced by a Hamiltonian vector field. In terms of the
Hamiltonian function $H(x,y)$, a point $(x,y)$ is conformal if two conditions
on the second partial derivatives hold:

$H_{xx}(x,y)=H_{yy}(x,y),\ H_{xy}(x,y)=0.$ (1)

There are two ways to interpret conditions (1). The first is to consider the
Hessian

$\mathit{Hess}(H)=\begin{pmatrix}H_{xx}&H_{yx}\\\ H_{xy}&H_{yy}\end{pmatrix}.$

This symmetric matrix has two orthogonal eigendirections, giving rise to two
fields of directions, and a conformal point is a singular point of these
fields, occurring when the Hessian is a scalar matrix.

This is similar to the fields of principal directions on a surface in
${\mathbb{R}}^{3}$. The singularities of theis field are the umbilic points of
the surface, the points where the two principal curvatures are equal. This
connection to classical differential geometry provides a context and it is one
of our motivations.

The famous Carathéodory conjecture asserts that a smooth closed surface in
${\mathbb{R}}^{3}$ homeomorphic to the sphere has at least two distinct
umbilic points. The sum of indices of the singular points of the field of
principal direction equals 2, the Euler characteristic of $S^{2}$. Since the
indices of singular points of line fields are half-integers, one concludes
that the algebraic number of umbilic points equals 4. See Figure 1 for the
umbilic points on a triaxial ellipsoid.

Figure 1: The four umbilic points, two of which are visible, and the lines of
curvature on an ellipsoid. This description is due to G. Monge, see [17] for a
historical overview.

The Carathéodory conjecture has a long and complicated history, see, e.g., [6,
17] and the references therein. The main approach to it in the literature is
local: conjecturally, the index of a singular point of the field of principal
directions does not exceed 1; this conjecture implies the Carathéodory
conjecture. Examples of umbilics with every index not exceeding 1 are known.

Another way to interpret conditions (1) is to consider the vector field

$V=(H_{yy}-H_{xx})\ \partial_{x}-2H_{xy}\ \partial_{y}$ (2)

whose zeros are the conformal points. The components of the vector field $V$
are the real and imaginary parts of $(\partial_{x}+i\partial_{y})^{2}H(x,y)$,
where $\partial_{x}+i\partial_{y}$ is the Cauchy-Riemann operator.

Loewner’s conjecture concerns the indices of isolated zeros of planar vector
fields $V_{n}$ whose two components are the real and imaginary parts of
$(\partial_{x}+i\partial_{y})^{n}H(x,y)$, where $H(x,y)$ is a smooth function.
We call these zeros Loewner points. Loewner’s conjecture states that the index
of an isolated zero of $V_{n}$ does not exceed $n$, see [19]. For $n=1$, this
is the well known statement that the index of an isolated zero of a gradient
vector field is not greater than 1.

In the case $n=2$, Loewner’s conjecture is intimately related with the
Carathéodory conjecture. In addition to differential geometry, Loewner’s
conjecture is related with hydrodynamics, see Section 3 of [15]. The
connection of conformal points to the Loewner conjecture is another motivation
for this study.

Let us formulate our results.

In Section 2 we prove that if $H$ is a smooth function in a simply connected
domain $D\subset{\mathbb{R}}^{2}$ such that $dH$ vanishes along the boundary
and the system (1) has no solutions on $\partial D$, then there are two
solutions inside $D$, counting with multiplicities. Under the same assumptions
on $H$, we extend this result to conformal points of the Hamiltonian vector
field of $H$ in a simply connected 2-dimensional domain with a general
Riemannian metric.

In Section 3 we extend this result to Loewner points: if a smooth function $H$
in a simply connected domain $D\subset{\mathbb{R}}^{2}$ is constant on the
boundary, its first $n-1$ normal derivatives vanish, and the respective vector
field $V_{n}$ has no zeros on $\partial D$, then there are $n$ Loevner points
inside $D$, counting with multiplicities.

In Section 4 we consider conformal points of symplectomorphisms of simply
connected domains $D\subset{\mathbb{R}}^{2}$. We prove that if a
symplectomorphism preserves $\partial D$ pointwise and is moderate then it
possess two conformal points, counted with multiplicity. The definition of a
moderate symplectomorphism is given in Section 4; informally speaking, a
moderate symplectomorphism is sufficiently close, with its first derivatives,
to the identity map.

In this sense, this result is akin to the first results on the existence of
symplectic fixed points, before the advent of the pseudo-holomorphic curves
and the Floer homology. We wonder whether the existence of conformal points
can be tackled using the techniques of the modern symplectic topology.

Section 5 concerns various interpretations of the Carathéodory conjecture and
its relation to conformal points. We show that, in a proper sense, umbilic
points are conformal points of a Hamiltonian vector field on $S^{2}$, and we
conjecture that an area preserving diffeomorphism of $S^{2}$ possesses at
least two distinct conformal points.

##  2 Conformal points of Hamiltonian vector fields

Let $D\subset{\mathbb{R}}^{2}$ be a connected and simply connected domain $D$
with non-empty interior, bounded by a smooth closed curve $\gamma$. Let
$H:D\to{\mathbb{R}}$ be a smooth function defined in a neighborhood of $D$ and
which is constant on $\gamma$. Thus, its Hamiltonian vector field
$X_{H}=H_{y}\partial_{x}-H_{x}\partial_{y}$ is tangent to the boundary
$\gamma$. We ask if the time-$\varepsilon$-flow $F_{\varepsilon}:D\to D$ of
$X_{H}$ has conformal points. That this is not true without further
assumptions is illustrated by the following example suggested to us by D.
Panov.

######  Example 2.1.

Let $\gamma$ be the ellipse $x^{2}+cy^{2}=1$ and $H(x,y)=x^{2}+cy^{2}$. If
$c\neq 1$, the Hamiltonian vector field $H_{y}\partial_{x}-H_{x}\partial_{y}$
has no conformal points, i.e., the infinitesimal symplectomorphism given by

$F:(x,y)\mapsto(x+\varepsilon H_{y},y-\varepsilon H_{x})=(x+2\varepsilon
cy,y-2\varepsilon x).$

has no conformal points. By continuity, it follows that, for sufficiently
small $\varepsilon\neq 0$, the time-$\varepsilon$-flow of this Hamiltonian
vector field is free from conformal points as well.

This example points towards a difference between this question and say a
fixed-point problem. Indeed, the time-$\varepsilon$-flow $F_{\varepsilon}$ of
an autonomous vector field $X$ is fixed-point free for all sufficiently small
$\varepsilon>0$ if and only if $X$ has no zeros. However, it is unclear to us
if the existence of a conformal point for the Hamiltonian vector field forces
the existence of conformal point for the time-$\varepsilon$-flow.
Nevertheless, the infinitesimal question is interesting on its own right and
links to the Carathéodory and Loewner conjectures as we mentioned in the
introduction. We return to the question for time-$\varepsilon$-flows in
Section 4.

######  Problem 2.2.

Is there an example of a Hamilton function $H:D\to{\mathbb{R}}$ such that the
Hamiltonian vector field has a conformal point but for all sufficiently small
$\varepsilon\neq 0$ the time-$\varepsilon$-flow of $H$ has no conformal point?

Let us now consider a general infinitesimal symplectomorphism

$F:(x,y)\mapsto(x+\varepsilon H_{y},y-\varepsilon H_{x})$

generated by $H:D\to{\mathbb{R}}$ with $H$ being constant on $\gamma$. The
Jacobian of $F$ is

$\begin{pmatrix}1+\varepsilon H_{xy}&\varepsilon H_{yy}\\\ -\varepsilon
H_{xx}&1-\varepsilon H_{xy}\end{pmatrix}.$

$F$ has conformal points, by definition, if this Jacobian matrix is an
infinitesimal similarity, that is, if and only if

$H_{xx}=H_{yy},\ H_{xy}=0,$

(equation (1), i.e., if the Hessian $\mathit{Hess}(H)$ of $H$ is a scalar
matrix.

Let $\gamma(t)$ be an arc length parameterization of the boundary of $D$, and
let $k(t)$ be the curvature of $\gamma(t)$. Assume that $\gamma$ is oriented
as the boundary of $D$ and let $n(t)$ be the inward unit normal vector field
along $\gamma$. A sufficiently small neighborhood of $\gamma$ is parameterized
as $\gamma(t)+sn(t),s\in(-{\varepsilon},{\varepsilon})$. Near $\gamma$ we
write the Hamiltonian function as $H(t,s)$.

######  Theorem 1.

Let the function $H:D\to{\mathbb{R}}$ be constant on the boundary $\gamma$ of
$D$. We assume that the system of equations (1) has isolated solutions in $D$,
none of which lie on the boundary of $D$. If $H$ satisfies one of the
following conditions

  1. i)

$H_{ss}(t,0)+k(t)H_{s}(t,0)$ does not vanish on $\gamma$,

  2. ii)

$H$ has zero normal derivative along $\gamma$,

then (1) has two solutions in the interior of $D$, counting with
multiplicities.

Before giving a proof, we make two comments. First, since $H$ is constant on
$\gamma$, one can write $H(t,s)=c+sh(t,s)$, where $h$ is a smooth function.
Then along the boundary $\gamma$ of $D$

$H_{ss}(t,0)+k(t)H_{s}(t,0)=2h_{s}(t,0)+k(t)h(t,0).$

If $\gamma$ has positive curvature, that is, $k>0$, one can always add a
constant to $h$ to ensure that $H_{ss}(t,0)+k(t)H_{s}(t,0)>0$ for all $t$.

Second, since $H$ is constant on $\gamma$, having zero normal derivative is
equivalent to $dH=0$ on $\gamma$, that is, to the Hamiltonian vector field
$H_{y}\partial_{x}-H_{x}\partial_{y}$ being identically zero on the boundary.

###### Proof.

We consider the vector field

$V=(H_{yy}-H_{xx})\ \partial_{x}-2H_{xy}\ \partial_{y}$

(equation (2). By assumption, $V$ has finitely many zeros in $D$, none of
which lie on the boundary. We want to show that $V$ has at least two zeros
inside $D$, counted with multiplicities. For that purpose we compute the
vector field $V$ along $\gamma$ in $(t,s)$-coordinates.

Let $\alpha(t)$ be the direction of $\gamma(t)$, that is,
$\gamma_{t}=(\cos\alpha,\sin\alpha)$. Then $k=\alpha_{t}$ and
$n=(-\sin\alpha,\cos\alpha)$. Combining the assumption that $H$ is constant on
$\gamma$, i.e., $H_{t}(t,0)=H_{tt}(t,0)=0$, with a calculation using the chain
rule yields

$\begin{split}V(t,0)=\;\;\;\;&[(H_{ss}(t,0)+k(t)H_{s}(t,0))\cos
2\alpha+2H_{st}(t,0)\sin 2\alpha]\ \partial_{x}\\\\[2.15277pt]
+&[(H_{ss}(t,0)+k(t)H_{s}(t,0))\sin 2\alpha-2H_{st}(t,0)\cos 2\alpha]\
\partial_{y},\end{split}$

or, in complex notation,

$V=(H_{ss}+kH_{s}-2iH_{st})e^{2i\alpha}.$

Under the first assumption, i.e., if $H_{ss}(t,0)+k(t)H_{s}(t,0)\neq 0$ for
all $t$, one has that $V$ has non-zero inner product with the vector field
$\cos(2\alpha)\partial_{x}+\sin(2\alpha)\partial_{y}$. The latter has rotation
number 2 along $\gamma$, hence the rotation number of $V$ also equals 2. The
Poincaré-Hopf theorem implies that $V$ has two zeros inside $D$, counting
multiplicities.

Under the second assumption, i.e., if $H$ has zero normal derivative along
$\gamma$, then $H_{s}(t,0)=H_{st}(t,0)=0$, and $V=H_{ss}e^{2i\alpha}$. By
assumption, $V$ has no zeros on the boundary, i.e., in this case, $H_{ss}$
does not vanish on the boundary. In particular, $V$ has again rotation number
2 along $\gamma$ and we conclude as above finishing the proof of Theorem 1.

Let us present an alternative argument in the second case. Instead of $V$ we
consider the Hessian matrix $\mathit{Hess}(H)$. At a point where this matrix
is not scalar, i.e., at a non-conformal point, it has two distinct eigenvalues
with orthogonal eigendirections. One may distinguish between the two
eigendirections by the value of the eigenvalues. Thus one obtains two
orthogonal line fields whose common singularities are the conformal points.

The rotation number of a line field along a closed curve is its number of the
complete turns divided by 2. In particular, the indices of singular points of
line fields are half-integers. This convention is consistent with the rotation
number for vector fields: when a non-oriented line makes a full turn, the
respective oriented line makes only half a turn. See Figure 2 for a
singularity of a line field with index $\tfrac{1}{2}$.

Figure 2: A singularity of a line field with index $\tfrac{1}{2}$.

The Poincaré-Hopf theorem holds also for line fields: the rotation number on
the boundary of $D$ equals the sum of indices of its singularities inside $D$.

As explained above, the second assumption is equivalent to $dH=0$ on the
boundary curve $\gamma$, i.e., $\gamma$ is a critical curve of the function
$H$. Therefore the tangent line to $\gamma$ is in the kernel of
$\mathit{Hess}(H)$, hence it is an eigendirection. If there are no conformal
points on the boundary, then $\gamma$ is an integral curve of one of the two
fields of eigendirections of $\mathit{Hess}(H)$. It follows that the rotation
number of this field around $\gamma$ is 1, and by the Poincaré-Hopf theorem it
has to have singular points inside $D$. Since their indices are half-integers,
one has at least two such points, counted with multiplicities. ∎

######  Remark 2.3.

We saw that the vector field $V$ rotates twice as fast as the field of
eigendirections of $\mathit{Hess}(H)$. This fact is known in the literature on
the Carathéodory conjecture, see, e.g., [15].

The next example shows that, under the assumptions of Theorem 1, one may have
exactly one conformal point (necessarily of multiplicity 2). We also discuss a
limitation of the above argument if we allow for conformal points on the
boundary of $D$.

######  Example 2.4.

Let $D$ be the closed unit disc, and let $H(\alpha,r)=r^{a}-\frac{a}{2}r^{2}$,
where $(\alpha,r)$ are polar coordinates on $D$ and $a>2$. This function is at
least of class $C^{2}$ (and if $a$ is an even integer, it is a polynomial), it
is constant and has zero normal derivative on the boundary $r=1$.

A calculation similar to the one in the proof of Theorem 1 yields

$V=\left(H_{rr}-\frac{1}{r}H_{r}\right)(\cos 2\alpha,\sin
2\alpha)=a(a-2)r^{a-2}(\cos 2\alpha,\sin 2\alpha).$

This vector field $V$ is non-singular in the punctured disc, and its only zero
is the origin.

We now consider $H$ as a function on
$D^{\prime}:=\\{z\in{\mathbb{C}}\mid|z-1|\leq 1\\}$, the unit disk around
$1\in{\mathbb{C}}$. Then $H$ has no conformal points in the interior but one
of multiplicity 2 on the boundary. Moreover, the vector field $V$ has
precisely one zero which lies on the boundary.

As in the alternative argument above, one is tempted to consider the singular
line field generated by the vector field $V$ (with one singularity in
$0\in\partial D^{\prime}$) since it has a continuous extension to a non-
singular line field along the boundary, namely the line field generated by the
vector field $(\cos 2\alpha,\sin 2\alpha)$. Now the conclusion from the
Poincaré-Hopf theorem seems to be that our example necessarily has two
conformal points in the interior of $D^{\prime}$ which is wrong as explained
above. The error lies in the fact that, while the line field indeed has a non-
singular continuous extension to the boundary given by $(\cos 2\alpha,\sin
2\alpha)$, there is no continuous extension to all of $D^{\prime}$.

It is unclear to us what the correct statement of Theorem 1 is if we allow for
zeros of $V$, i.e., conformal points, on the boundary.

######  Problem 2.5.

Is there a Hamiltonian function $H:D\to{\mathbb{R}}$ such that the Hamiltonian
vector field has only one conformal point with multiplicity 1 on the boundary
of $D$?

The circle of questions makes sense on other surfaces. Since there are no
topological obstructions on the torus coming from the Poincaré-Hopf theorem we
pose the following problem.

######  Problem 2.6.

Let $H(x,y)$ be a function on the torus ${\mathbb{R}}^{2}/L$ where $L$ is a
lattice. Does the respective Hamiltonian vector field necessarily have
conformal points? What about Hamiltonian symplectomorphisms?

We point out that elements from $SL(2,{\mathbb{Z}})$ acts on
${\mathbb{R}}^{2}/{\mathbb{Z}}^{2}$ by area-preserving maps and, except for
the $\pm$identity, none of them have conformal points. However, again except
for the identity, none of them are Hamiltonian diffeomorphisms, i.e.,
generated by a Hamiltonian function on ${\mathbb{R}}^{2}/{\mathbb{Z}}^{2}$.

Finally, we extend case ii) of Theorem 1 to domains $D\subset{\mathbb{R}}^{2}$
as above but equipped with a general Riemannian metric. This metric induces a
symplectic form and a conformal structure, and we may again consider conformal
points of Hamiltonian vector fields.

######  Theorem 2.

Let $dH=0$ on the boundary $\gamma$ so that the Hamiltonian vector field
$X_{H}$ of $H$ has no conformal points on $\gamma$. Then $X_{H}$ has two
conformal points in the interior of the domain, counting with multiplicities.

###### Proof.

Let the metric be given in the conformal form
$\tfrac{1}{g(x,y)}(dx^{2}+dy^{2})$, where $g$ is a positive function. Then the
symplectic form is $\tfrac{1}{g(x,y)}dx\wedge dy$, and therefore the
Hamiltonian vector field of a function $H(x,y)$ is
$g(x,y)(H_{y}\partial_{x}-H_{x}\partial_{y})$. A point $(x,y)$ is conformal
for the Hamiltonian vector field if and only if

$(gH_{x})_{x}=(gH_{y})_{y},\ (gH_{y})_{x}=-(gH_{x})_{y}.$

As before, consider the vector field

$V=[(gH_{y})_{y}-(gH_{x})_{x}]\ \partial_{x}-[(gH_{x})_{y}+(gH_{y})_{x}]\
\partial_{y}.$

If $dH=0$ on the boundary then, on $\gamma$, one has

$V=g[(H_{yy}-H_{xx})\ \partial_{x}-2H_{xy}\ \partial_{y}].$

This field still has rotation number 2 along $\gamma$, and the rest of the
proof is the same as that of Theorem 1. ∎

##  3 Loewner points

Consider the vector fields $V_{n}$ whose components are the real and imaginary
parts of the function $(\partial_{x}+i\partial_{y})^{n}H(x,y)$. One has
$V_{1}=\nabla H$, and $-V_{2}$ is the vector field $V$ that we considered
above. Recall that zeros of $V_{n}$ are termed Loewner points.

We extend case (ii) of Theorem 1 to Loewner points. Let $D$ be as in Theorem
1.

######  Theorem 3.

Assume that the function $H$ is constant on the boundary $\gamma$ and has
vanishing $j$-th normal derivatives along $\gamma$ for $j=1,\ldots,n-1$.
Moreover, assume that the vector field $V_{n}$ has only isolated zeros in $D$
and no zeros on the boundary. Then $V_{n}$ has $n$ zeros in the interior of
$D$, counting with multiplicities.

###### Proof.

We use the same coordinates $(t,s)$ as before, and we again calculate the
vector field $V_{n}$ along $\gamma$. The chain rule implies that

$\partial_{x}+i\partial_{y}=e^{i\alpha(t)}\left(\frac{1}{1-sk(t)}\
\partial_{t}+i\partial_{s}\right)=:e^{i\alpha(t)}{\mathcal{D}}$

where we recall that $\alpha_{t}=k$. Moreover, we observe

$\displaystyle f$
$\displaystyle=\mathcal{D}(e^{in\alpha(t)}f)-e^{in\alpha(t)}\mathcal{D}f$
$\displaystyle=(\mathcal{D}e^{in\alpha(t)})f+e^{in\alpha(t)}\mathcal{D}f-e^{in\alpha(t)}\mathcal{D}f$
$\displaystyle=ne^{in\alpha(t)}\frac{i\alpha_{t}(t)}{1-sk(t)}f,$

that is, the first-order differential operator $\mathcal{D}$ satisfies the
commutation relation

$[{\mathcal{D}},e^{in\alpha(t)}]=ne^{in\alpha(t)}\frac{ik(t)}{1-sk(t)}=:ne^{in\alpha(t)}F(t,s),$

where the right hand side is the zero-order operator of multiplying by a
function. Lemma 3.1 below then implies that

$(\partial_{x}+i\partial_{y})^{n}=e^{ni\alpha}({\mathcal{D}}+(n-1)F)({\mathcal{D}}+(n-2)F)\ldots({\mathcal{D}}+F){\mathcal{D}}.$

Using the specific form of $\mathcal{D}$ one has

$({\mathcal{D}}+(n-1)F)({\mathcal{D}}+(n-2)F)\ldots({\mathcal{D}}+F){\mathcal{D}}H(t,s)=\sum_{1\leq
j+k\leq n}c_{jk}(t,s)\partial_{t}^{j}\partial_{s}^{k}H(t,s),$

where $c_{jk}(t,s)$ are some functions. Evaluating the right hand side on the
boundary $\gamma$ and using the assumption that the first $n-1$ normal
derivatives of $H$ vanish on $\gamma$, i.e., $\partial_{s}^{k}H(t,0)=0$ for
$j=1,\ldots,n-1$, the only term that survives is
$i^{n}(\partial_{s}^{n}H)(t,0)$. It follows that on the boundary one has

$V_{n}=i^{n}(\partial_{s}^{n}H)(t,0)\ e^{in\alpha(t)},$

and the proof concludes the same way as that of Theorem 1. ∎

######  Lemma 3.1.

Let $\mathcal{D}$ be a first order differential operator on real-valued
functions in several variables. We assume that $\mathcal{D}$ satisfies the
commutation relation $[\mathcal{D},a^{n}]=na^{n}b$, $n\in\mathbb{N}$, for some
functions $a$ and $b$, where $a^{n}$ resp. $na^{n}b$ are understood as
multiplication operators. Then

$(a\mathcal{D})^{n}=a^{n}\prod_{k=1}^{n}(\mathcal{D}+(n-k)b)$

holds for all $n\in\mathbb{N}$.

###### Proof.

We proceed by induction. The claim is clearly true for $n=1$. For the
induction step let $f$ be some function. To simplify the following computation
we set $g:=\prod_{k=1}^{n}(\mathcal{D}+(n-k)b)f$. Then

$\displaystyle(a\mathcal{D})^{n+1}f$
$\displaystyle=(a\mathcal{D})(a\mathcal{D})^{n}f$
$\displaystyle=a\mathcal{D}\left(a^{n}\prod_{k=1}^{n}(\mathcal{D}+(n-k)b)f\right)$
$\displaystyle=a\mathcal{D}(a^{n}g)$
$\displaystyle=a\big{(}a^{n}\mathcal{D}g+na^{n}bg\big{)}$
$\displaystyle=a^{n+1}(\mathcal{D}+nb)g$
$\displaystyle=a^{n+1}(\mathcal{D}+nb)\prod_{k=1}^{n}(\mathcal{D}+(n-k)b)f$
$\displaystyle=a^{n+1}\prod_{k=0}^{n}(\mathcal{D}+(n-k)b)f$
$\displaystyle=a^{n+1}\prod_{k=1}^{n+1}(\mathcal{D}+(n+1-k)b)f$

completes the proof by induction. We used the induction assumption in the
second equation and the commutation relation in the fourth equation. ∎

##  4 Conformal points of moderate symplectomorphisms

We now consider an extension of the second case of Theorem 1 to
symplectomorphisms of a domain in ${\mathbb{R}}^{2}$ satisfying a certain
condition to be defined below.

We continue to assume that $D$ is as above, that is, a connected and simply
connected domain $D$, bounded by a smooth closed curve $\gamma$. Consider a
symplectomorphism $F:D\to D$, and let $G\subset D\times
D\subset{\mathbb{C}}^{2}$ be its graph. We equip ${\mathbb{C}}^{2}$ with the
symplectic form $\omega_{1}\ominus\omega_{2}$, where $\omega_{1}$ and
$\omega_{2}$ are the pullbacks of the standard area form in
${\mathbb{C}}={\mathbb{R}}^{2}$ to the first and the second summands. Let
${\mathbb{C}}\cong\Delta\subset{\mathbb{C}}^{2}$ be the diagonal, the graph of
the identity map. Graphs of symplectomorphisms are Lagrangian submanifolds.

Let $L:(z_{1},z_{2})\mapsto(w_{1},w_{2})$ be the complex linear map
${\mathbb{C}}^{2}\to{\mathbb{C}}^{2}$ given by the formulas

$w_{1}=\frac{z_{1}+z_{2}}{2},\ w_{2}=i(z_{1}-z_{2}).$

We denote $z_{j}=x_{j}+iy_{j},w_{j}=u_{j}+iv_{j}$, $j=1,2$.

######  Lemma 4.1.

The map $L$ is a symplectomorphism between
$({\mathbb{C}}^{2},\omega_{1}\ominus\omega_{2})$ and
$(T^{*}\Delta\cong{\mathbb{C}}^{2},d\lambda)$ where
$\lambda=u_{2}du_{1}+v_{2}dv_{1}$.

###### Proof.

The required equality

$dx_{1}\wedge dy_{1}-dx_{2}\wedge dy_{2}=du_{2}\wedge du_{1}+dv_{2}\wedge
dv_{1}$

is verified by a direct calculation. ∎

Consider a symplectomorphism $F$ whose graph $G\subset{\mathbb{C}}^{2}$ is
transverse to the fibers of the cotangent bundle $T^{*}\Delta$, that is, $G$
needs to be transverse to the projection $\mathrm{pr}:{\mathbb{C}}^{2}\cong
T^{*}\Delta\to\Delta$ which in the identification of Lemma 4.1 simply is the
map

$(z_{1},z_{2})\mapsto\mathrm{pr}(z_{1},z_{2})=\left(\frac{z_{1}+z_{2}}{2},\frac{z_{1}+z_{2}}{2}\right).$

We call such symplectomorphisms moderate.

In particular, $F$ is moderate if it is $C^{1}$-close to the identity. More
precisely, one has the next result. Let the symplectomorphism of $D$ be given
by $(x,y)\mapsto F(x,y)=(f(x,y),g(x,y))$, in particular,
$f_{x}g_{y}-f_{y}g_{x}=1$.

######  Lemma 4.2.

The graph $G$ of the symplectomorphism $F$ is transverse to the fibers of
$T^{*}\Delta$ near a point $(z,F(z))\in G$ if and only if
$2+f_{x}(z)+g_{y}(z)\neq 0$ with $z=x+iy$. In particular, $F$ is moderate if
and only if $2+f_{x}+g_{y}\neq 0$ for all $(x,y)\in D$.

######  Remark 4.3.

In other words the graph $G$ a symplectomorphism $F$ is transverse to the
fibers of $T^{*}\Delta$ if its Jacobian in $\mathrm{SL(2,{\mathbb{R}})}$ is
never negative parabolic. It easy to see that, e.g., the graph of a shear
$(x,y)\mapsto(-x\pm y,-y)$ indeed is not transverse to the fibers.

###### Proof.

One has $G=\\{(x,y,f,g)\mid(x,y)\in D\\}$, and its tangent space is spanned by
the vectors $(1,0,f_{x},g_{x})$ and $(0,1,f_{y},g_{y})$. The kernel of the
projection on the diagonal is spanned by the vectors $(1,0,-1,0)$ and
$(0,1,0,-1)$. These two spaces are not transverse if and only if

$\det\begin{pmatrix}1&0&f_{x}&g_{x}\\\ 0&1&f_{y}&g_{y}\\\ 1&0&-1&0\\\
0&1&0&-1\end{pmatrix}=0.$

This determinant equals $(1+f_{x})(1+g_{y})-f_{y}g_{x}=2+f_{x}+g_{y}$, and the
result follows. ∎

Let $\bar{D}\subset\Delta$ be the diagonal embedding of $D$, i.e., the graph
of $\mathrm{id}_{D}$, and denote by $\bar{\gamma}$ its boundary.

######  Lemma 4.4.

Assume that $F$ is moderate and is the identity on the boundary $\gamma$. Then
$\mathrm{pr}:G\to\bar{D}$ is a diffeomorphism.

###### Proof.

That $\mathrm{pr}$ is a mop onto $\bar{D}$ is a fairly general topological
argument and does not use that $G$ is transverse to the fibers.

Indeed, assume that $B\subset{\mathbb{R}}^{n}$ is smooth $n$-dimensional
manifold with boundary homeomorphic to the closed unit ball and
$\Phi:B\to{\mathbb{R}}^{n}$ is a continuous map with $\Phi|_{\partial
B}=\mathrm{id}$. We claim that $B\subset\Phi(B)$ holds.

Assume for a contradiction that we find a point $p\in B\setminus\Phi(B)$. Then
we may consider the map $\varphi:\partial B\hookrightarrow
B\stackrel{{\scriptstyle\Phi}}{{\to}}{\mathbb{R}}^{n}\setminus\\{p\\}\stackrel{{\scriptstyle
r}}{{\simeq}}\partial B$ where the map
$r:{\mathbb{R}}^{n}\setminus\\{p\\}\to\partial B$ is a retraction. Here we use
that $p\in B$, in fact, $p\in\mathring{B}$ since $\Phi|_{\partial
B}=\mathrm{id}$. Moreover, $\Phi|_{\gamma}$ is not null-homotopic again since
$\Phi|_{\partial B}=\mathrm{id}$. On the other hand $\Phi$ factors through the
contractible set $B$, i.e., is null-homotopic, and thus so is
$\Phi|_{\gamma}$. This contradiction shows that $B\subset\Phi(B)$.

Applying this to $\Phi:B:=\bar{D}\to{\mathbb{C}}$ with
$\Phi(z,z):=\mathrm{pr}(z,F(z))$ we conclude
$\bar{D}\subset\mathrm{pr}(G)=\Phi(\bar{D})$.

Now assume that $F$ is in addition transverse to $\mathrm{pr}$. Then
$\mathrm{pr}$ is a local diffeomorphism. Let us prove that
$\mathrm{pr}(G)\subset\bar{D}$. This would imply that $\mathrm{pr}$ is a
covering, but since $\bar{D}$ is simply connected, it must be a
diffeomorphism.

Choose a point $(w,w)\in\mathrm{pr}(G)\setminus\bar{D}\subset\Delta$ having
maximal distance to some fixed point $o\in\bar{D}$. Choose $(z,F(z))\in G$
with $\mathrm{pr}(z,F(z))=(w,w)$. Since $\mathrm{pr}$ is a local
diffeomorphism the point $(z,F(z))$ needs to be on the boundary of $G$:
otherwise we could move its projection $(w,w)=\mathrm{pr}(z,F(z))$ further
away from $o$. Since $G$ is diffeomorphic to $D$ and $F$ is the identity on
$\partial D$ we actually conclude $(z,F(z))=(z,z)=(w,w)\in\partial\bar{D}$.

In conclusion, the points of $\mathrm{pr}(G)$ of maximal distance to $o$ are
on the boundary of $\bar{D}$, in particular, $\mathrm{pr}(G)\subset\bar{D}$,
as claimed. ∎

######  Remark 4.5.

We point out that in the previous proof of $\mathrm{pr}(G)\subset\bar{D}$ we
did neither use that $F$ is a symplectomorphism nor that $D$ is simply
connected. Moreover, this statement has the following geometric
interpretation. The mid-point map $D\ni
z\mapsto\tfrac{1}{2}(z+F(z))\in{\mathbb{C}}$ takes values inside $D$ if $F$ is
moderate and the identity on the boundary. For non-convex domains $D$ this is
geometrically not directly clear to us.

We are ready to extend Theorem 1 to moderate symplectomorphisms.

######  Theorem 4.

Let $F$ be a moderate symplectomorphism of $D$ which is the identity on the
boundary $\gamma$. Assume that $F$ has no conformal points on $\gamma$. Then
$F$ possesses at least two conformal points inside $D$, counted with
multiplicity.

###### Proof.

A point $(x,y)$ is conformal if and only if

$f_{x}(x,y)=g_{y}(x,y),\ f_{y}(x,y)=-g_{x}(x,y).$

We claim that, counted with multiplicity, $f_{y}+g_{x}$ and $g_{y}-f_{x}$
vanish simultaneously at at least two points. We shall reduce this to case ii)
in Theorem 1.

In the notation of Lemma 4.1, $x=x_{1},y=y_{1}$ and, accordingly, we write $u$
for $u_{1}$ and $v$ for $v_{1}$. Thus $(u,v)$ are coordinates in $\bar{D}$,
and Lemma 4.4 implies that $(x,y)\mapsto(u,v)$ is a diffeomorphism. Therefore
we may consider $f_{x}(x,y)-g_{y}(x,y)$ and $f_{y}(x,y)+g_{x}(x,y)$ as
functions of $u$ and $v$. We “pack” these two functions in a vector filed

$V=(f_{x}(x,y)-g_{y}(x,y))\ \partial u+(f_{y}(x,y)+g_{x}(x,y))\ \partial v.$

Our goal is to show that $V$ has two zeros inside $\bar{D}$. For this purpose,
we calculate $V$ on the boundary $\bar{\gamma}$.

The graph $G$ of $F$ is Lagrangian and, by assumption, transverse to the
fibers of $T^{*}\Delta$. Thus the graph is a section of the cotangent bundle
given by the differential of a function $H:\bar{D}\to{\mathbb{R}}$. Since
$F=(f,g)$ is the identity on $\gamma$, one has $dH=0$ on $\bar{\gamma}$.

Using the coordinate change from Lemma 4.1 we obtain the following formulas
involving the graph $G=\\{(x,y,f,g)\mid(x,y)\in D\\}$ and $dH$

$\displaystyle 2x$ $\displaystyle=2u+H_{v},\quad$ $\displaystyle 2f$
$\displaystyle=2u-H_{v},$ $\displaystyle 2y$ $\displaystyle=2v-H_{u},\quad$
$\displaystyle 2g$ $\displaystyle=2v+H_{u}.$

Taking differentials, it follows that

$\displaystyle 2dx$ $\displaystyle=2du+H_{uv}du+H_{vv}dv,\quad$ $\displaystyle
2f_{x}dx+2f_{y}dy$ $\displaystyle=2du-H_{uv}du-H_{vv}dv,$ $\displaystyle 2dy$
$\displaystyle=2dv-H_{uu}du-H_{uv}dv,\quad$ $\displaystyle 2g_{x}dx+2g_{y}dy$
$\displaystyle=2dv+H_{uu}du+H_{uv}dv.$

Substitute the two left equations into the right ones and equate the resulting
1-forms to obtain

$\displaystyle f_{x}\big{(}4-H_{uv}^{2}+H_{uu}H_{vv}\big{)}$
$\displaystyle=4+H_{uv}^{2}-H_{uu}H_{vv}-4H_{uv},$ $\displaystyle
f_{y}\big{(}4-H_{uv}^{2}+H_{uu}H_{vv}\big{)}$ $\displaystyle=-4H_{vv},$
$\displaystyle g_{x}\big{(}4-H_{uv}^{2}+H_{uu}H_{vv}\big{)}$
$\displaystyle=4H_{uu},$ $\displaystyle
g_{y}\big{(}4-H_{uv}^{2}+H_{uu}H_{vv}\big{)}$
$\displaystyle=4+H_{uv}^{2}-H_{uu}H_{vv}+4H_{uv}.$

Denote by $t$ the variable along $\bar{\gamma}$. Then deriving $H_{u}=0$ resp.
$H_{v}=0$ along $\bar{\gamma}$ we obtain from the chain rule

$\displaystyle H_{uu}u_{t}+H_{vu}v_{t}$ $\displaystyle=0,\
H_{uv}u_{t}+H_{vv}v_{t}$ $\displaystyle=0.$

It follows that $H_{uu}H_{vv}=H_{uv}^{2}$, and we obtain

$\displaystyle f_{y}+g_{x}=H_{uu}-H_{vv},\quad g_{y}-f_{x}=2H_{uv}$

along $\bar{\gamma}$. Therefore, along $\bar{\gamma}$, one has

$V=(H_{uu}-H_{vv})\ \partial u+2H_{uv}\ \partial v,$

which takes us to the second case of Theorem 1. Its proof shows exactly the
desired assertion. ∎

##  5 Carathéodory conjecture and complex points of Lagrangian surfaces

###  5.1 Geometry of the space of oriented lines

In this section we survey the geometry of the space of oriented lines in
${\mathbb{R}}^{3}$ and the space of oriented non-parameterized geodesics in
$S^{3}$ and $H^{3}$. We refer to [1, 8, 9, 10] and Section 5.6 of [14] for the
material of this and the subsequent sections.

To start with, the space ${\mathcal{G}}(M)$ of oriented non-parameterized
geodesics of a Riemannian manifold $M$ carries a symplectic structure that is
obtained from the canonical symplectic structure of $T^{*}M$ by symplectic
reduction. This is true locally, and if the space of geodesics is a smooth
manifold, then this manifold is symplectic.

If $N\subset M$ is a cooriented hypersurface, then the normal geodesics to $N$
provide a Lagrangian immersion of $N$ into ${\mathcal{G}}(M)$. Denote the
image of this immersion by ${\mathcal{L}}(N)$.

In the Euclidean case, one has ${\mathcal{G}}({\mathbb{R}}^{n})=T^{*}S^{n-1}$,
in particular, ${\mathcal{G}}({\mathbb{R}}^{3})=T^{*}S^{2}$. For hyperbolic
3-space, the space of geodesics is also symplectomorphic to $T^{*}S^{2}$. If
$M=S^{3}$, then the geodesics are the great circles, and
${\mathcal{G}}(S^{3})=G_{+}(2,4)$, the Grassmannian of oriented 2-dimensional
subspaces in $R^{4}$. One has: $G_{+}(2,4)=S^{2}\times S^{2}$, see, e.g., [7].

Given a line $\ell\subset{\mathbb{R}}^{3}$, its train
$S(\ell)\subset{\mathcal{G}}({\mathbb{R}}^{3})$ consists of the lines that
intersect $\ell$. The train is a singular hypersurface. In the tangent space
$T_{\ell}{\mathcal{G}}({\mathbb{R}}^{3})$, one obtains a quadratic cone. A
similar construction works for $H^{3}$ and $S^{3}$.

Next, ${\mathcal{G}}({\mathbb{R}}^{3})$ carries a complex structure. Let
$\ell$ be an oriented line in ${\mathbb{R}}^{3}$. The $90^{\circ}$ rotation of
the ambient space about $\ell$ induces an automorphism of the tangent space
$T_{\ell}{\mathcal{G}}({\mathbb{R}}^{3})$; one obtains an almost complex
structure, and this structure is integrable. Similarly for $H^{3}$ and
$S^{3}$.

Note that this complex structure is not compatible with the symplectic
structure. For example, the lines through a fixed point form a Lagrangian
sphere which is also a complex curve.

Furthermore, ${\mathcal{G}}({\mathbb{R}}^{3})$ carries a Kähler structure
whose metric has signature $(+,+,-,-)$. Given a surface
$N\subset{\mathbb{R}}^{3}$, the restriction of this metric to
${\mathcal{L}}(N)$ is either Lorentz or null. The latter happens at the
complex points of ${\mathcal{L}}(N)$.

The intersection of the quadratic cone in $T{\mathcal{G}}({\mathbb{R}}^{3})$
with $T{\mathcal{L}}(N)$ is the light cone of the restriction of the Kähler
metric to ${\mathcal{L}}(N)$. The two light directions on ${\mathcal{L}}(N)$
correspond to the two principal directions on the surface $N$. Likewise for
$H^{3}$ and $S^{3}$.

###  5.2 Umbilics as conformal points

The differential equation of the principal directions is given by

$\beta_{xy}dx^{2}+(\beta_{yy}-\beta_{xx})dxdy-\beta_{xy}dy^{2}=0,$

where $\beta(x,y)$ is the Bonnet function, defined in a special Bonnet chart
on the surface, see [5]. Representing the principal directions as the
eigendirections of a quadratic form is similar to our alternative proof of
Theorem 1 using the Hessian matrix $\mathit{Hess}(H)$.

Let $H:S^{2}\to{\mathbb{R}}$ be a smooth function. Consider the respective
Hamiltonian vector field on the sphere.

######  Lemma 5.1.

A point $s\in S^{2}$ is conformal for the Hamiltonian vector field of $H$ if
and only if the second jet $j^{2}H(s)$ equals the second jet of a first
spherical harmonic, that is, the restriction to $S^{2}$ of an affine function
$d+ax+by+cz$ in ${\mathbb{R}}^{3}$.

###### Proof.

A point is conformal for an area preserving diffeomorphism if its derivative
sends circles to congruent circles, that is, it is an isometry. An orientation
preserving isometry of $S^{2}$ is a rotation about an axis, and its
Hamiltonian is the restriction to $S^{2}$ of an affine function $d+ax+by+cz$
(the axis being parallel to the vector $(a,b,c)$). ∎

Recall that the support function of a convex surface
$N\subset{\mathbb{R}}^{3}$ is a function $H:S^{2}\to{\mathbb{R}}$ that equals
the (signed) distance $H(s)$ from the origin to the tangent plane of $N$ whose
oriented normal is $s\in S^{2}$. The differential $dH$ defines a Lagrangian
section of $T^{*}S^{2}$; using its identification with
${\mathcal{G}}({\mathbb{R}}^{3})$, we obtain the Lagrangian surface
${\mathcal{L}}(N)$.

######  Lemma 5.2.

The umbilic points of $N$ correspond to points $s\in S^{2}$ where $j^{2}H(s)$
is the second jet of a first spherical harmonic.

###### Proof.

The umbilic points are the points where a sphere is 2nd order tangent to $N$,
that is, where $j^{2}H$ equals the second jet of the support function of a
sphere. But the support functions of the spheres are the first harmonics (with
$d>0$ being the radius of the sphere). ∎

Lemmas 5.1 and 5.2 imply that the umbilic points of a convex surface $N$ are
the conformal points of the Hamiltonian vector field of the support function
of $N$. We are led to the following generalization of the Carathéodory
conjecture.

######  Conjecture 1.

An area preserving diffeomorphism of $S^{2}$ possesses at least two distinct
conformal points.

######  Example 5.3.

Consider the particular case of Theorem 2 for the spherical metric. Assume
that the support function has vanishing differential along a curve
$\gamma\subset S^{2}$. Then the respective surface $N$ is tangent to a sphere
along $\gamma$. The normal lines to $N$ along $\gamma$ form a cone, a
developable surface.

It is a classical fact that the normals to a surface along a curve form a
developable surface if and only if the curve is a line of curvature, see,
e.g., [18]. (Indeed, the normal lines at infinitesimally close points of a
line of curvature intersect, and the locus of these intersection points is a
curve in ${\mathbb{R}}^{3}$. Then the surface comprising the normals is the
tangent developable of this space curve, a developable surface.)

It follows that $\gamma$ is a closed line of curvature of $N$, and hence there
there is an umbilic point in each of the two components of its complement.

Let us also interpret umbilics in terms of the space of oriented lines. This
interpretation works equally well in the spherical and hyperbolic geometries.
Let $N\subset{\mathbb{R}}^{3}$ be a cooriented surface.

######  Lemma 5.4.

The umbilic points of $N$ correspond to the complex points of
${\mathcal{L}}(N)\subset{\mathcal{G}}({\mathbb{R}}^{3})$.

###### Proof.

A point of $N$ is umbilic if and only if the quadratic surface that
approximates $N$ at this point up to the second derivatives is a sphere, say
$S$. That is, ${\mathcal{L}}(N)$ and ${\mathcal{L}}(S)$ are tangent at the
respective point. But ${\mathcal{L}}(S)$ is a complex curve in
${\mathcal{G}}({\mathbb{R}}^{3})$, and the result follows. ∎

Thus (a slightly generalized) formulation of the Carathéodory conjecture is
that a Lagrangian section of $T^{*}S^{2}$ possesses at least two distinct
complex points.

######  Remark 5.5.

There exists a substantial literature concerning complex points of real
surfaces. For example, consider an oriented closed surface $M$ embedded in
${\mathbb{C}}^{2}$. The generic complex points of $M$ are classified into four
types: positive or negative (the orientation of the tangent plane is complex
or not), and elliptic or hyperbolic, see [4, 3] and [2], pp. 137–139. One has
the topological relations:

$e_{+}-h_{+}=e_{-}-h_{-}=\frac{1}{2}\chi(M),$

where $\chi(M)$ is the Euler characteristic. In particular, if
$M\subset{\mathbb{C}}^{2}$ is a sphere in general position, then there exist
at least two elliptic complex points.

Acknowledgements. Many thanks to D. Panov and A. Petrunin for useful
discussions.

PA acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through Germany’s Excellence Strategy EXC-2181/1 -
390900948 (the Heidelberg STRUCTURES Excellence Cluster), the Transregional
Colloborative Research Center CRC/TRR 191 (281071066). ST was supported by NSF
grant DMS-2005444 and by a Mercator fellowship within the CRC/TRR 191, and he
thanks the Heidelberg University for its hospitality.

## References

  * [1] D. Alekseevsky, B. Guilfoyle, W. Klingenberg. On the geometry of spaces of oriented geodesics. Ann. Global Anal. Geom. 40 (2011), 389–409.
  * [2] D. Bennequin. Entrelacements et équations de Pfaff. Astérisque, 107–108, Soc. Math. France, Paris, 1983.
  * [3] E. Bishop. Differentiable manifolds in complex Euclidean space. Duke Math. J. 32 (1965), 1–21.
  * [4] S.-S. Chern, E. Spanier. A theorem on orientable surfaces in four-dimensional space. Comment. Math. Helv. 25 (1951), 205–209.
  * [5] G. Darboux. Sur la forme des lignes de courbure dans la voisinage d’un ombilic. Note 7, Leçons sur la Théorie générale des Surfaces IV. Gauthier-Villars, Paris, 1896.
  * [6] M. Ghomi, R. Howard. Normal curvatures of asymptotically constant graphs and Carathéodory’s conjecture. Proc. Amer. Math. Soc. 140 (2012), 4323–4335.
  * [7] H. Gluck, F. Warner. Great circle fibrations of the three-sphere. Duke Math. J. 50 (1983), 107–132.
  * [8] B. Guilfoyle, W. Klingenberg. On the space of oriented affine lines in ${\mathbb{R}}^{3}$. Arch. Math. 82 (2004), 81–84.
  * [9] B. Guilfoyle, W. Klingenberg. Generalised surfaces in ${\mathbb{R}}^{3}$. Math. Proc. R. Ir. Acad. 104A (2004), 199–209.
  * [10] B. Guilfoyle, W. Klingenberg. An indefinite Kähler metric on the space of oriented lines. J. London Math. Soc. 72 (2005), 497–509.
  * [11] C. Gutierrez, J. Sotomayor. Lines of curvature, umbilic points and Carathéodory conjecture. Resenhas 3 (1998), 291–322.
  * [12] H.F. Lai. Characteristic classes of real manifolds immersed in complex manifolds. Trans. Amer. Math. Soc. 172 (1972), 1–33.
  * [13] C. Loewner. A topological characterization of a class of integral operators. Ann. of Math. 49 (1948), 316–332
  * [14] V. Ovsienko, S. Tabachnikov. Projective differential geometry old and new. From the Schwarzian derivative to the cohomology of diffeomorphism groups.  Cambridge Univ. Press, Cambridge, 2005.
  * [15] B. Smyth, F. Xavier. Real solvability of the equation $\partial_{\bar{z}}^{2}=\rho g$ and the topology of isolated umbilics. J. Geom. Anal. 8 (1998), 655–671.
  * [16] B. Smyth, F. Xavier. Eigenvalue estimates and the index of Hessian fields. Bull. London Math. Soc. 33 (2001), 109–112.
  * [17] J. Sotomayor, R. Garcia. Lines of curvature on surfaces, historical comments and recent developments. São Paulo J. Math. Sci. 2 (2008), 99–143.
  * [18] D. Struik. Lectures on Classical Differential Geometry. Addison-Wesley Press, Inc., Cambridge, Mass., 1950.
  * [19] C. Titus. A proof of a conjecture of Loewner and of the conjecture of Caratheodory on umbilic points. Acta Math. 131 (1973), 43–77.
  * [20] https://mathoverflow.net/questions/354451.