\section*{Introduction}

The theory of optimal transport \cite{Vil03,Vil09} has experienced rapid growth in recent years with applications in diverse fields across pure and applied mathematics. Along with this growth came a lot of interest in extending the methods of optimal transport beyond the scope of its original formulation as an optimization problem for the transport cost between two probability measures.

One such extension deals with ``quantum spaces'', where the probability measures are replaced by density matrices or density operators. Most of the work dealing with quantum optimal transport in this sense can be grouped into one of the following two categories. The first approach (see e.g. \cite{CM14a,CM17a,MM17,CGT18,DR19,BV20,CM20}) relies on a noncommutative analog of the Benamou--Brenier formulation \cite{BB00} of the Wasserstein distance for probability measures on Euclidean space
\begin{equation*}
W_2^2(\mu,\nu)=\inf\left\lbrace\int_0^1\int_{\mathbb{R}^n}\abs{v_t}^2\,d\rho_t\,dt: \rho_0=\mu,\rho_1=\nu,\dot\rho_t+\nabla\cdot(\rho_t v_t)=0\right\rbrace.
\end{equation*}
This approach has proven fruitful in applications to noncommutative functional inequalities, similar in spirit to the heuristics known as Otto calculus \cite{CM17a,CM20,DR20,WZ20}.

The second approach (see e.g. \cite{NGT15,GMP16,PCVS19,Duv20,PMTL20,PT21}) seeks to find a suitable noncommutative analog of the Monge--Kantorovich formulation \cite{Kan42} of the Wasserstein distance via couplings (or transport plans):
\begin{equation*}
W_p^p(\mu,\nu)=\inf\left\lbrace\int_{X\times X}d^p(x,y)\,d\pi(x,y):(\mathrm{pr}_1)_\#\pi=\mu,(\mathrm{pr}_2)_\#\pi=\nu\right\rbrace.
\end{equation*}
This approach also allows to consider a quantum version of Monge--Kantorovich problem for arbitrary cost functions. So far, possible connections between these two approaches in the quantum world stay elusive.

The focus of this article lies on the noncommutative transport distance $\W$ introduced in the first approach. More precisely, we prove a dual formula that is a noncommutative analog of the expression of the classical $L^2$-Wasserstein distance in terms of subsolutions of the Hamilton--Jacobi equation \cite{OV00,BGL01}
\begin{equation*}
W_2^2(\mu,\nu)=\frac 1 2 \inf\left\lbrace\int_{\mathbb{R}^n} u_1\,d\mu-\int_{\mathbb{R}^n}u_0\,d\nu: \dot u_t+\frac 1 2\abs{\nabla u_t}^2\leq 0\right\rbrace.
\end{equation*}
This result yields a noncommutative version of the dual formula obtained independently by Erbar, Maas and the author \cite{EMW19} and Gangbo, Li and Mou \cite{GLM19} for the Wasserstein-like transport distance on graphs. In fact, we prove a dual formula that is not only valid for the metric $\W$, but also for the entropic regularization recently introduced by Becker--Li \cite{BL21}.

With the notation introduced in the next section, the main result of this article reads as follows.
\begin{theorem*}
Let $\sigma\in M_n(\IC)$ be an invertible density matrix and $(P_t)$ an ergodic quantum Markov semigroup on $M_n(\IC)$ that satisfies the $\sigma$-DBC. The entropic regularization $\W_\epsilon$ of noncommutative transport distance induced by $(P_t)$ satisfies the following dual formula:
\begin{equation*}
\frac 1 2\W_\epsilon^2(\rho_0,\rho_1)=\sup\{\tau(A(1)\rho_1-A(0)\rho_0)\mid A\in \HJB^1_{\epsilon}\}.
\end{equation*}
\end{theorem*}
Here $\HJB^1_{\epsilon}$ stands for the set of all Hamilton--Jacobi--Bellmann subsolutions, a suitable noncommutative variant of solutions of the differential inequality
\begin{equation*}
\dot u(t)+\frac 1 2\abs{\nabla u(t)}^2-\epsilon \Delta u(t)\leq 0.
\end{equation*}

Other metrics similar to $\W$ also occur in the literature, most notably the one called the ``anticommutator case'' in \cite{CGT18,CGGT20,BL21}. In \cite{Wir18,CM20}, a class of such metrics was studied in a systematic way, and our main theorem applies in fact to this wider class of metrics. For the anticommutator case, this duality formula was obtained before in \cite{CGGT20}.

There are still some very natural questions left open. For one, we do not discuss the existence of optimizers. While for the primal problem this follows from a standard compactness argument, this question is more delicate for the dual problem, even when dealing with probability densities on discrete spaces instead of density matrices, and one has to relax the problem to obtain maximizers (see \cite[Sections 6--7]{GLM19}).

Another interesting direction would be to extend the duality result from matrix algebras to infinite-dimensional systems. While a definition of the metric $\W$ for quantum Markov semigroups on semi-finite von Neumann algebras is available \cite{Hor18,Wir18}, the problem of duality seems to be much harder to address. Even for abstract diffusion semigroups, the best known result only shows that the primal distance is the upper length distance associated with the dual distance and leaves the question of equality open \cite[Proposition 10.11]{AES16}.

\subsection*{Acknowledgments}
The author wants to thank Jan Maas for helpful comments. He acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 716117) and from the Austrian Science Fund (FWF) through project F65.

\section{Setting and basic definitions}

In this section we introduce basic facts and definitions about quantum Markov semigroups that will be used later on. In particular, we review the definition of the noncommutative transport distance from \cite{CM17a} and its entropic regularization introduced in \cite{BL21}. Our notation mostly follows \cite{CM17a,CM20}.

Let $M_n(\IC)$ denote the complex $n\times n$ matrices and let $\A$ be a unital $\ast$-subalgebra of $M_n(\IC)$. Let $\A_h$ denote the self-adjoint part of $\A$. We write $\tau$ for the normalized trace on $M_n(\IC)$ and $\H_\A$ for the Hilbert space formed by equipping $\A$ with the GNS inner product
\begin{equation*}
\langle\cdot,\cdot\rangle_{\H_A}\colon \A\times\A\to\IC,\,(A,B)\mapsto \tau(A^\ast B).
\end{equation*}
The adjoint of a linear operator $\mathscr{K}\colon \H_\A\to\H_\A$ is denoted by $\mathscr{K}^\dagger$.

We write $\S(\A)$ for the set of all density matrices on $\A$, that is, all positive elements $\rho\in\A$ with $\tau(\rho)=1$. The subset of invertible density matrices is denoted by $\S_+(\A)$.

A quantum Markov semigroup (QMS) on $\A$ is a family $(P_t)_{t\geq 0}$ of linear operators on $\A$ that satisfy the following conditions:
\begin{itemize}
\item $P_t$ is unital and completely positive for every $t\geq 0$,
\item $P_0=\mathrm{id}_\A$, $P_{s+t}=P_s P_t$ for all $s,t\geq 0$,
\item $t\mapsto P_t$ is continuous.
\end{itemize}
We consider a QMS $(P_t)$ on $\A$ which extends to a QMS on $M_n(\IC)$ satisfying the $\sigma$-DBC for some density matrix $\sigma\in \S_+(A)$, that is,
\begin{equation*}
\tau((P_t A)^\ast B \sigma)=\tau(A^\ast (P_t B)\sigma)
\end{equation*}
for $A,B\in\A$. Let $\L$ denote the generator of $(P_t)$, that is, the linear operator on $\A$ given by
\begin{equation*}
\L(A)=\lim_{t\searrow 0}\frac{P_t A-A}{t}.
\end{equation*
We further assume that $(P_t)$ is \emph{ergodic} (or primitive), that is, the kernel of $\L$ is one-dimensional.

By Alicki's theorem \cite[Theorem 3]{Ali76}, \cite[Theorem 3.1]{CM17a} there exists a finite set $\J$, real numbers $\omega_j$ for $j\in \J$ and $V_j\in M_n(\IC)$ for $j\in \J$ with the following properties:
\begin{itemize}
\item $\tau(V_j^\ast V_k)=\delta_{jk}$ for $j,k\in\J$,
\item $\tau(V_j)=0$ for $j\in\J$,
\item $\{V_j\mid j\in\J\}=\{V_j^\ast\mid j\in\J\}$,
\item $\sigma V_j \sigma^{-1}=e^{-\omega_j}V_j$ for $j\in \J$
\end{itemize}
such that
\begin{equation*}
\L(A)=\sum_{j\in \J}\left(e^{-\omega_j/2}V_j^\ast[A,V_j]-e^{\omega_j/2}[A,V_j]V_j^\ast\right)
\end{equation*}
for $A\in \A$.

The numbers $\omega_j$ are called Bohr frequencies of $\L$ and are uniquely determined by $(P_t)$.

The matrices $V_j$ are not uniquely determined by $(P_t)$ and $\sigma$, but in the following we will fix a set $\{V_j\mid j\in\J\}$ that satisfies the preceding conditions.

Let
\begin{equation*}
\H_{\A,\J}=\bigoplus_{j\in \J}\H_\A^{(j)},
\end{equation*}
where $\H_\A^{(j)}$ is a copy of $\H_\A$ for $j\in\J$. We write $\partial_j$ for $[V_j,\cdot\,]$ and
\begin{equation*}
\nabla\colon \H_\A\to\H_{\A,\J},\,\nabla(A)=(\partial_j(A))_{j\in\J}.
\end{equation*}
We write $\div$ for the adjoint of $-\nabla$, that is,
\begin{equation*}
\div=-\sum_{j\in\J}\partial_j^\dagger.
\end{equation*}

For $X\in\A_+$ and $\alpha\in \IR$ define
\begin{equation*}
[X]_\alpha\colon \H_\A\to \H_\A,\,[X]_\alpha(A)=\int_0^1 e^{\alpha(s-1/2)}X^s A X^{1-s}\,ds.
\end{equation*}
Given $\vec\alpha=(\alpha_j)_{j\in\J}$, we define
\begin{equation*}
[X]_{\vec\alpha}\colon \H_{\A,\J}\to\H_{\A,\J},\,(A_j)_{j\in\J}\mapsto ([X]_{\alpha_j}A_j)_{j\in\J}.
\end{equation*}

Now let $\vec\omega=(\omega_j)_{j\in\J}$ with the Bohr frequencies $\omega_j$ of $\L$. For $\epsilon\geq 0$ we write $\CE_\epsilon(\rho_0,\rho_1)$ for the set of all pairs $(\rho,\bV)$ such that $\rho\in H^1([0,1];\S_+(\A))$ with $\rho(0)=\rho_0$, $\rho(1)=\rho_1$, $\bV\in L^2([0,1];\H_{\A,\J})$ and
\begin{equation*}
\dot \rho(t) +\div\bV(t)=\epsilon \L^\dagger\rho(t)
\end{equation*}
for a.e. $t\in [0,1]$.

We define a metric $\W_\epsilon$ on $\S_+(\A)$ by
\begin{equation*}
\W_\epsilon^2(\rho_0,\rho_1)=\inf_{(\rho,\bV)\in \CE_\epsilon(\rho_0,\rho_1)}\int_0^1 \langle \bV(t),[\rho(t)]_{\vec\omega}^{-1}\bV(t)\rangle\,dt.
\end{equation*}
A standard mollification argument shows that the infimum can equivalently be taken over $(\rho,\bV)\in \CE_\epsilon(\rho_0,\rho_1)$ with $\rho\in C^\infty([0,1];\S_+(\A))$.

For $\epsilon=0$, this is the noncommutative transport distance $\W$ introduced in \cite{CM17a} (as distance function associated with a Riemannian metric on $\S(\A)_+$), and for $\epsilon>0$, this is the entropic regularization of $\W$ introduced in \cite{BL21}.

By a substitution one can reformulate the minimization problem for $\W_\epsilon$ in such a way that the constraint becomes independent from $\epsilon$. For that purpose define the relative entropy of $\rho\in \S_+(\A)$ with respect to $\sigma$ by
\begin{equation*}
D(\rho\|\sigma)=\tau(\rho(\log\rho-\log \sigma))
\end{equation*}
and the Fisher information of $\rho\in \S_+(\A)$ by
\begin{equation*}
\I(\rho)=\langle [\rho]_{\vec\omega}\nabla (\log \rho-\log\sigma),\nabla(\log\rho-\log\sigma)\rangle_{\H_{\A,\J}}.
\end{equation*}

According to \cite[Theorem 1]{BL21}, one has
\begin{align*}
\W_\epsilon^2(\rho_0,\rho_1)&=\inf_{(\rho,\bW)\in \CE_0(\rho_0,\rho_1)}\int_0^1 (\langle \bW(t),[\rho(t)]_{\vec\omega}^{-1}\bW(t)\rangle+\epsilon^2\I(\rho(t)))\,dt\\
&\quad+2\epsilon(D(\rho_1\|\sigma)-D(\rho_0\|\sigma)).
\end{align*}

The metric $\W$ is intimately connected to the relative entropy and therefore well-suited to study its decay properties along the quantum Markov semigroup. For other applications, variants of the metric $\W$ have also proven useful (e.g. \cite{CGT18,CGGT20}), for which the operator $[\rho]_{\vec\omega}$ is replaced. A systematic framework of these metrics has been developed in \cite{Wir18,CM20}. It can be conveniently phrased in terms of so-called operator connections.

Let $H$ be an infinite-dimensional Hilbert space. A map $\Lambda\colon B(H)_+\times B(H)_+\to B(H)_+$ is called an \emph{operator connection} \cite{KA80} if
\begin{itemize}
\item $A\leq C$ and $B\leq D$ imply $\Lambda(A,B)\leq \Lambda(C,D)$,
\item $C\Lambda(A,B)C\leq \Lambda(CAC,CBC)$,
\item $A_n\searrow A$, $B_n\searrow B$ imply $\Lambda(A_n,B_n)\searrow \Lambda(A,B)$.
\end{itemize}
For example, for every $\alpha\in \IR$ the map
\begin{equation*}
\Lambda_\alpha\colon (A,B)\mapsto \int_0^1 e^{\alpha(s-1/2)}A^s B^{1-s}\,ds
\end{equation*}
is an operator connection.

It can be shown that every operator connection $\Lambda$ satisfies
\begin{equation*}
U^\ast \Lambda(A,B)U=\Lambda(U^\ast AU,U^\ast BU)
\end{equation*}
for $A,B\in B(H)_+$ and unitary $U\in B(H)$ \cite[Section 2]{KA80}. Embedding $\IC^n$ into $H$, one can view $A,B\in M_n(\IC)$ as bounded linear operators on $H$, and the unitary invariance of $\Lambda$ ensures that $\Lambda(A,B)$ does not depend on the embedding of $\IC^n$ into $H$.

For $X\in \A$ define
\begin{align*}
L(X)&\colon \H_{\A}\to \H_{\A},\,A\mapsto XA\\
R(X)&\colon \H_{\A}\to \H_{\A},\,A\mapsto AX.
\end{align*}
With this notation we can write
\begin{equation*}
[\rho]_\Lambda=\Lambda(L(\rho),R(\rho)).
\end{equation*}

Since $L(X)$ and $R(X)$ commute, we have
\begin{equation}\label{eq:spec_rep}
\Lambda(L(X),R(X))A=\sum_{k,l=1}^n \Lambda(\lambda_k,\lambda_l)E_k A E_l\tag{$\ast$}
\end{equation}
for $X\in \A_+$ and $A\in \H_\A$, where $(\lambda_k)$ are the eigenvalues of $X$ and $E_k$ the corresponding spectral projections.

More generally let $\vec\Lambda=(\Lambda_j)_{j\in \J}$ be a family of operator connections and define
\begin{align*}
[\rho]_{\Lambda_j}&=\Lambda_j(L(\rho),R(\rho)),\\
[\rho]_{\vec \Lambda}&=\bigoplus_{j\in \J}[\rho]_{\Lambda_j}.
\end{align*}
Then one can define a distance $\W_{\vec \Lambda}$ by
\begin{align*}
\W_{\vec\Lambda,\epsilon}(\rho_0, \rho_0)^2=\inf_{(\rho,\bV)\in\CE_\epsilon(\rho_0,\rho_1)}\int_0^1 \langle [\rho(t)]_{\vec\Lambda}^{-1}\bV(t),\bV(t)\rangle_{\H_{\A,\J}}\,dt.
\end{align*}
.

If $\Lambda_j=\Lambda_{\omega_j}$ as above, then we retain the original metric $\W_\epsilon$, while for $\Lambda_j(A,B)=\frac 1 2 (A+B)$ (and $\epsilon=0$) one obtains the distance studied in \cite{CGT18,CGGT20}.

Later we will make the additional assumption that $\Lambda_{j^\ast}(A,B)=\Lambda_j(B,A)$. It follows from the representation theorem of operator means \cite{KA80} that the class of metrics $\W_{\vec\Lambda,0}$ with $\vec\Lambda$ subject to this symmetry condition is exactly the class of metrics satisfying Assumptions 7.2 and 9.5 in \cite{CM20}.

For technical reasons, it can be useful to allow for curves of density matrices that are not necessarily invertible. For this purpose, we make the following convention: If $\mathcal{K}\colon \H_{\A,\J}\to\H_{\A,\J}$ is a positive operator and $\bV\in\H_{\A,\J}$, we define
\begin{equation*}
\langle \bV,\mathcal{K}^{-1}\bV\rangle_{\H_{\A,\J}}=\begin{cases}\langle\mathcal{K}\bW,\bW\rangle_{\H_{\A,\J}}&\text{if }\bV\in (\ker\mathcal{K})^\perp,\mathcal{K}\bW=\bV,\\\infty&\text{otherwise}.\end{cases}
\end{equation*}
Since $(\ker \mathcal{K})^\perp=\operatorname{ran}\mathcal{K}$ and $\mathcal{K}$ is injective on $(\ker\mathcal{K})^\perp$, the element $\bW$ in this definition exists and is unique. Moreover, this convention is clearly consistent with the usual definition if $\mathcal{K}$ is invertible.

Alternatively, as a direct consequence of the spectral theorem, $\langle \bV,\mathcal K^{-1}\bV\rangle_{\H_{\A,\J}}$ can be expressed as
\begin{equation*}
\langle \bV,\mathcal{K}^{-1}\bV\rangle_{\H_{\A,\J}}=\sum_{k=1}^m \frac 1{\lambda_k}\abs{\langle\bV,\bW_k\rangle_{\H_{\A,\J}}}^2,
\end{equation*}
where $\lambda_1,\dots,\lambda_m$ are the eigenvalues of $\mathcal{K}$ and $\W_1,\dots,W_m$ an orthonormal basis of corresponding eigenvectors.

\begin{lemma}\label{lem:monotone_conv}
If $\mathcal K_n\colon \H_{\A,\J}\to\H_{\A,\J}$, $n\in \mathbb N$, are positive invertible operators that converge monotonically decreasing to $\mathcal{K}$, then
\begin{equation*}
\langle \bV,\mathcal K_n^{-1}\bV\rangle_{\H_{\A,\J}}\nearrow \langle \langle \bV,\mathcal K^{-1}\bV\rangle_{\H_{\A,\J}}
\end{equation*}
for all $\bV\in\H_{\A,\J}$.
\end{lemma}
\begin{proof}
From the spectral expression it is easy to see that
\begin{equation*}
\langle\bV,\mathcal{K}_n^{-1}\bV\rangle_{\H_{\A,\J}}=\sup_{\delta>0}\langle \bV,(\mathcal{K}_n+\delta)^{-1}\bV\rangle_{\H_{\A,\J}}
\end{equation*}
and the same for $\mathcal{K}_n$ replaced by $\mathcal{K}$. Moreover, since $\mathcal{K}_n\searrow \mathcal{K}$, we have $(\mathcal{K}_n+\delta)^{-1}\nearrow (\mathcal{K}+\delta)^{-1}$. Thus
\begin{align*}
\langle \bV,\mathcal{K}^{-1}\bV\rangle_{\H_{\A,\J}}&=\sup_{\delta>0}\langle \bV,(\mathcal{K}+\delta)^{-1}\bV\rangle_{\H_{\A,\J}}\\
&=\sup_{\delta>0}\sup_{n\in\mathbb N}\langle\bV,(\mathcal{K}_n+\delta)^{-1}\bV\rangle_{\H_{\A,\J}}\\
&=\sup_{n\in\mathbb N}\sup_{\delta>0}\langle\bV,(\mathcal{K}_n+\delta)^{-1}\bV\rangle_{\H_{\A,\J}}\\
&=\sup_{n\in\mathbb N}\langle \bV,\mathcal{K}_n^{-1}\bV\rangle_{\H_{\A,\J}}.
\end{align*}
Since $(\langle\bV,\mathcal{K}_n^{-1}\bV\rangle_{\H_{\A,\J}})$ is monotonically increasing, this settles the claim.
\end{proof}

Write $\CE_{\epsilon}^\prime(\rho_0,\rho_1)$ for the set of all pairs $(\rho,\bV)$ such that $\rho\in H^1([0,1];\S(\A))$ with $\rho(0)=\rho_0$, $\rho(1)=\rho_1$, $\bV\in L^2([0,1];\H_{\A,\J})$ and
\begin{equation*}
\dot \rho(t)+\div\bV(t)=\epsilon\L^\dagger\rho(t)
\end{equation*}
for a.e. $t\in[0,1]$. The only difference to the definition of $\CE_\epsilon(\rho_0,\rho_1)$ is that $\rho(t)$ is not assumed to be invertible.

\begin{proposition}\label{prop:relaxed_paths_CE}
For $\rho_0,\rho_1\in \S_+(\A)$ we have
\begin{equation*}
\W_{\vec\Lambda,\epsilon}^2(\rho_0,\rho_1)=\inf_{(\rho,\bV)\in \CE_\epsilon^\prime(\rho_0,\rho_1)}\int_0^1 \langle \bV(t),[\rho(t)]_{\vec\Lambda}^{-1}\bV(t)\rangle\,dt.
\end{equation*}
\end{proposition}
\begin{proof}
It suffices to show that every curve $(\rho,\bV)\in \CE_\epsilon^\prime(\rho_0,\rho_1)$ can be approximated by curves in $\CE_\epsilon(\rho_0,\rho_1)$ such that the action integrals converge.

For that purpose let
\begin{equation*}
\rho^\delta\colon [0,1]\to \S_+(\A),\,t\mapsto\begin{cases}(1-t)\rho_0+t1_\A&\text{if }t\in[0,\delta],\\
(1-\delta)\rho((1-2\delta)^{-1}(t-\delta))+\delta 1_\A&\text{if }t\in(\delta,1-\delta),\\
t\rho_1+(1-t)1_\A&\text{if }t\in[1-\delta,1].
\end{cases}
\end{equation*}
Since $(P_t)$ is assumed to be ergodic, for $t\in[0,\delta]$ there exists $\bV^\delta(t)=\nabla U^\delta(t)$ such that
\begin{equation*}
\dot \rho^\delta(t)+\div\bV^\delta(t)=1-\rho_0+\div \bV^\delta(t)=\epsilon (1-t)\L^\dagger\rho_0=\epsilon \L^\dagger\rho^\delta(t)
\end{equation*}
and
\begin{equation*}
\norm{\bV^\delta(t)}_{\H_{\A,\J}}\leq C\norm{\epsilon(1-t)\L^\dagger\rho_0-1+\rho_0}_{\H_\A}
\end{equation*}
with a constant $C>0$ (depending only on the spectral gap of $\L$). In particular, it is bounded independent of $\delta$.

Moreover, if $\lambda$ is the smallest eigenvalue of $\rho_0$, which is strictly positive by assumption, then $\rho^\delta(t)\geq ((1-t)\lambda+t)1_\A\geq \lambda 1_\A$.

Thus
\begin{align*}
\int_0^\delta \langle \bV^\delta(t),[\rho^\delta(t)]_{\vec\Lambda}^{-1}\bV^\delta(t)\rangle_{\H_{\A,\J}}\,dt&\leq \int_0^\delta \langle \bV^\delta(t),[\lambda 1_\A]_{\vec\Lambda}^{-1}\bV^\delta(t)\rangle_{\H_{\A,\J}}\,dt\\
&\leq \norm{[\lambda 1_\A]_{\vec \Lambda}^{-1}}\int_0^\delta\norm{\bV^\delta(t)}_{\H_{\A,\J}}^2\,dt\\
&\to 0
\end{align*}
as $\delta\to 0$. Similarly one can show
\begin{equation*}
\lim_{\delta\to 0}\int_{1-\delta}^1 \langle \bV^\delta(t),[\rho^\delta(t)]_{\vec\Lambda}^{-1}\bV^\delta(t)\rangle_{\H_{\A,\J}}\,dt=0.
\end{equation*}
By the same argument as above, for a.e. $t\in (\delta,1-\delta)$ there exists a unique gradient $\bW^\delta(t)$ such that
\begin{equation*}
\div \bW^\delta(t)=-\frac{2\delta\epsilon}{1-2\delta}\L^\dagger \rho((1-2\delta)^{-1}(t-\delta))
\end{equation*}
and
\begin{equation*}
\norm{\bW^\delta(t)}_{\H_{\A,\J}}\leq \frac{2\delta\epsilon}{1-2\delta}\norm{\L^\dagger\rho((1-2\delta)^{-1}(t-\delta))}_{\H_{\A,\J}}.
\end{equation*}
Since $\rho\in H^1([0,1];\S(\A))\subset C([0,1];\S(\A))$, the norm on the right side is bounded independent of $\delta$, so that
\begin{equation*}
\norm{\bW^\delta(t)}_{\H_{\A,\J}}\leq \tilde C\delta
\end{equation*}
with a constant $\tilde C>0$ independent of $\delta$. As $\rho^\delta(t)\geq \delta 1_\A$ for $t\in (\delta,1-\delta)$, this implies
\begin{align*}
\int_\delta^{1-\delta}\langle\bW^\delta(t), [\rho^\delta(t)]_{\vec\Lambda}^{-1}\bW^\delta(t)\rangle_{\H_{\A,\J}}\,dt&\leq \frac 1 \delta\int_\delta^{1-\delta}\langle \bW^\delta(t),[1_\A]_{\vec\Lambda}^{-1}\bW^\delta(t)\rangle_{\H_{\A,\J}}\,dt\\
&\leq \tilde C \norm{[1_\A]_{\vec\Lambda}^{-1}}\delta\\
&\to 0
\end{align*}
as $\delta\to 0$.

With
\begin{equation*}
\bV^\delta(t)=\frac 1{1-2\delta}\bV((1-2\delta)^{-1}(t-\delta))+\bW^\delta(t)
\end{equation*}
we have
\begin{equation*}
\dot \rho^\delta(t)+\div \bV^\delta(t)=\epsilon \L\rho^\delta(t).
\end{equation*}
Furthermore,
\begin{align*}
&\quad\int_\delta^{1-\delta}\left\langle\bV\left(\frac{t-\delta}{1-2\delta}\right),[\rho^\delta(t)]_{\vec\Lambda}^{-1}\bV\left(\frac{t-\delta}{1-2\delta}\right)\right\rangle_{\H_{\A,\J}}\,dt\\
&=\frac 1{1-\delta}\int_\delta^{1-\delta}\left\langle\bV\left(\frac{t-\delta}{1-2\delta}\right),\left[\rho\left(\frac{t-\delta}{1-2\delta}\right)+\frac{\delta}{1-\delta}\right]_{\vec\Lambda}^{-1}\bV\left(\frac{t-\delta}{1-2\delta}\right)\right\rangle_{\H_{\A,\J}}\,dt\\
&=\frac{1-2\delta}{1-\delta}\int_0^1 \langle \bV(s),\left[\rho(s)+\frac{\delta}{1-\delta}\right]_{\vec\Lambda}^{-1}\bV(s)\rangle_{\H_{\A,\J}}\,ds,
\end{align*}
where we used the substitution $s=(1-2\delta)^{-1}(t-\delta)$ in the last step.

By Lemma \ref{lem:monotone_conv} and the monotone convergence theorem we obtain
\begin{align*}
&\quad\;\lim_{\delta\to 0}\int_\delta^{1-\delta}\left\langle\bV\left(\frac{t-\delta}{1-2\delta}\right),[\rho^\delta(t)]_{\vec\Lambda}^{-1}\bV\left(\frac{t-\delta}{1-2\delta}\right)\right\rangle_{\H_{\A,\J}}\,dt\\
&=\int_0^1 \langle \bV(s),[\rho(s)]_{\vec\Lambda}^{-1}\bV(s)\rangle_{\H_{\A,\J}}\,ds.
\end{align*}

Together with the convergence result for $\bW^\delta$ from above, this implies
\begin{equation*}
\int_\delta^{1-\delta}\langle \bV^\delta(t),[\rho^\delta(t)]_{\vec \Lambda}^{-1}\bV^\delta(t)\rangle_{\H_{\A,\J}}\,dt\to \int_0^1\langle \bV(t),[\rho(t)]_{\vec\Lambda}^{-1}\bV(t)\rangle_{\H_{\A,\J}}\,dt.
\end{equation*}
Altogether we have shown
\begin{equation*}
\lim_{\delta\to 0}\int_0^{1}\langle \bV^\delta(t),[\rho^\delta(t)]_{\vec \Lambda}^{-1}\bV^\delta(t)\rangle_{\H_{\A,\J}}\,dt= \int_0^1\langle \bV(t),[\rho(t)]_{\vec\Lambda}^{-1}\bV(t)\rangle_{\H_{\A,\J}}\,dt.\qedhere
\end{equation*}
\end{proof}

\section{Real subspaces}

Since the proof of the main result relies on convex analysis methods for real Banach spaces, we need to identify suitable real subspaces for our purposes. For $\A$ this is simply $\A_h$, but for $\H_{\A,\J}$ this is less obvious and will be done in the following.

For $j\in\J$ denote by $j^\ast$ the unique index in $\J$ such that $V_j^\ast=V_{j^\ast}$. Let $\tilde \H_ {\A}^{(j)}$ be the linear span of $\{X\partial_j A\mid A,X\in\A\}$, and define a linear map $J\colon \tilde \H_{\A}^{(j)}\to \tilde \H_{\A}^{(j^\ast)}$ by
\begin{align*}
J(X\partial_j A)=\partial_{j^\ast} (A^\ast)X^\ast.
\end{align*}
By the product rule, $(\partial_j A)X$ also belongs to $\tilde \H_{\A,\J}^{(j)}$ and $J((\partial_j A)X)=X^\ast \partial_{j^\ast}(A^\ast)$.

\begin{lemma}
The map $J$ is anti-unitary.
\end{lemma}
\begin{proof}
For $A,B,X,Y\in\A$ we have
\begin{align*}
\langle J(X\partial_j A),J(Y\partial_j B)\rangle_{\H_\A}&=\tau(X(AV_j-V_A)(V_j^\ast B^\ast-B^\ast V_j^\ast) Y^\ast)\\
&=\tau((B^\ast V_j^\ast-V_j^\ast B^\ast)Y^\ast X(V_j A-AV_j))\\
&=\langle Y\partial_j B,X\partial_j A\rangle.\qedhere
\end{align*}
\end{proof}

Let
\begin{equation*}
\Hh_{\A,\J}=\{\bV\in \bigoplus_{j\in\J}\tilde\H_{\A}^{(j)}\mid J(\bV_j)=\bV_{j^\ast}\}.
\end{equation*}
By the previous lemma, $\Hh_{\J,\A}$ is a real Hilbert space.

\begin{lemma}\label{lem:grad_real}
Let $(\Lambda_j)_{j\in\J}$ be a family of operator connections such that $\Lambda_{j^\ast}(A,B)=\Lambda_j(B,A)$ for all $j\in\J$. If $A\in \A_h$ and $\rho\in \S(\A)$, then $\nabla A,[\rho]_{\vec\Lambda}\nabla A\in \Hh_{\A,\J}$.
\end{lemma}
\begin{proof}
For $\nabla A$ the statement follows directly from the definitions. For $[\rho]_{\vec\Lambda} \nabla A$ first note that
\begin{equation*}
J\Lambda(L(\rho),R(\rho))=\Lambda(R(\rho),L(\rho))J
\end{equation*}
as a consequence of the spectral representation (\ref{eq:spec_rep}).

Thus
\begin{align*}
J([\rho]_{\Lambda_j}\partial_j A)&=J\Lambda_j(L(\rho),R(\rho))\partial_j A\\
&=\Lambda_j(R(\rho),L(\rho))J\partial_j A\\
&=\Lambda_{j^\ast}(L(\rho),R(\rho))\partial_{j^\ast}A.\qedhere
\end{align*}
\end{proof}

\section{Duality}

In this section we prove the duality theorem announced in the introduction. Our strategy follows the same lines as the proof in the commutative case in \cite{EMW19}. It crucially relies on the Rockefellar--Fenchel duality theorem quoted below. Throughout this section we fix a quantum Markov semigroup with generator $\L$ satisfying the $\sigma$-DBC for some $\sigma\in \S_+(\A)$ and a family $(\Lambda_j)_{j\in\J}$ of operator connections such that $\Lambda_{j^\ast}(A,B)=\Lambda_j(B,A)$ for all $j\in \J$.

We need the following definition for the constraint of the dual problem. Here and in the following we write
\begin{equation*}
\langle \bV,\bW\rangle_\rho=\langle \bV,[\rho]_{\vec\Lambda}\bW\rangle_{\H_{\A,\J}}
\end{equation*}
for $\bV,\bW\in \H_{\A,\J}$ and $\rho\in \A_+$.

\begin{definition}\label{def:HJB}
A function $A\in H^1((0,T);\A_h)$ is said to be a \emph{Hamilton--Jacobi--Bellmann subsolution} if for a.e. $t\in (0,T)$ we have
\begin{align*}
\tau((\dot A(t)+\epsilon \L A(t)) \rho)+\frac 1 2 \norm{\nabla A(t)}_{\rho}^2\leq 0\qquad\text{for all }\rho\in \S(A).
\end{align*}
The set of all Hamilton--Jacobi--Bellmann subsolutions is denoted by $\HJB_{\vec\Lambda,\epsilon}$.
\end{definition}

Our proof will establish equality between the primal and dual problem, but before we begin, let us show that one inequality is actually quite easy to obtain.

\begin{proposition}
For all $\rho_0,\rho_1\in\S_+(\A)$ we have
\begin{equation*}
\sup_{A\in\HJB_{\vec\Lambda,\epsilon}}\tau(A(1)\rho_1-A(1)\rho_0)\leq\frac 1 2 \inf_{(\rho,\bV)\in\CE_\epsilon(\rho_0,\rho_1)}\int_0^1 \langle \bV(t),[\rho(t)]_{\vec\Lambda}^{-1}\bV(t)\rangle_{\H_{\A,\J}}\,dt.
\end{equation*}
\end{proposition}
\begin{proof}
For $A\in \HJB_{\vec\Lambda,\epsilon}$ and $(\rho,\bV)\in \CE_\epsilon(\rho_0,\rho_1)$ we have
\begin{align*}
\tau(A(1)\rho_1-A(0)\rho_0)&=\int_0^1 \tau(\dot A(t) \rho(t)+A(t)\dot\rho(t))\,dt\\
&\leq-\int_0^1 (\epsilon\tau((\L A(t))\rho(t))+\frac 1 2\norm{\nabla A(t)}_{\rho(t)}^2)\,dt\\
&\quad+\int_0^1(\langle \nabla A(t),\bV(t)\rangle_{\H_{\A,\J}}+\epsilon\tau(A(t)\L^\dagger\rho(t)))\,dt\\
&= \int_0^1\langle  [\rho(t)]_{\vec\omega}^{1/2}\nabla A(t),[\rho(t)]_{\vec\Lambda}^{-1/2}\bV(t)\rangle_{\H_{\A,\J}}\,dt\\
&\quad-\frac 1 2 \int_0^1\langle [\rho(t)]_{\vec\Lambda}^{1/2}\nabla A(t),[\rho(t)]_{\vec \Lambda}^{1/2}\nabla A(t)\rangle_{\H_{\A,\J}}\,dt\\
&\leq\frac 1 2\int_0^1\langle\bV(t),[\rho(t)]_{\vec\Lambda}^{-1}\bV(t)\rangle_{\H_{\A,\J}}\,dt,
\end{align*}
where we used $A\in \HJB_{\vec\Lambda,\epsilon}$ and $(\rho,\bV)\in \CE_\epsilon(\rho_0,\rho_1)$ for the first inequality and Young's inequality for the second inequality.
\end{proof}

To prove actual equality, our crucial tool is the Rockefellar--Fenchel duality theorem (see e.g. \cite[Theorem 1.9]{Vil03}, which we quote here for the convenience of the reader. Recall that if $E$ is a (real) normed space, the Legendre--Fenchel transform $F^\ast$ of a proper convex function $F\colon E\to \IR\cup\{\infty\}$ is defined by
\begin{equation*}
F^\ast\colon E^\ast\to\mathbb{R}\cup\{\infty\},\,F^\ast(x^\ast)=\sup_{x\in E}(\langle x^\ast,x\rangle-F(x)).
\end{equation*}

\begin{theorem}\label{thm:Rockefellar-Fenchel}
Let $E$ be a real normed space and $F,G\colon E\to \IR\cup\{\infty\}$ proper convex functions with Legendre--Fenchel transforms $F^\ast,G^\ast$. If there exists $z_0\in E$ such $G$ is continuous at $z_0$ and $F(z_0),G(z_0)<\infty$, then
\begin{equation*}
\sup_{z\in E}(-F(z)-G(z))=\min_{z^\ast\in E^\ast}(F^\ast(z^\ast)+G^\ast(-z^\ast)).
\end{equation*}
\end{theorem}

Before we state the main result, we still need the following useful inequality.

\begin{lemma}\label{lem:deriv_conv_fct}
For any operator connection $\Lambda$ the map
\begin{equation*}
f_{\Lambda}\colon \A_{++}\to B(\H_A),\,A\mapsto [A]_\Lambda
\end{equation*}
is smooth and its Fréchet derivative satisfies
\begin{align*}
d f_\Lambda(B) A\geq f_\Lambda(A)
\end{align*}
for $A,B\in \A_{++}$ with equality if $A=B$.
\end{lemma}
\begin{proof}
Smoothness of $f_\Lambda$ is a consequence of the representation theorem of operator connections [Theorem 3.4]\cite{KA80}. For the claim about the Fréchet derivative first note that $f_\Lambda$ is concave \cite[Theorem 3.5]{KA80}. Therefore $d^2 f_\Lambda(X)[Y,Y]\leq 0$ for all $X\in \A_{++}$ and $Y\in \A_h$ by \cite[Proposition 2.2]{Han97}.

The fundamental theorem of calculus implies
\begin{align*}
(d f_\Lambda(A)-d f_\Lambda(B))(A-B)&=\int_0^1 d^2 f_\Lambda(tA+(1-t)B)[A-B,A-B]\,dt\\
&\leq 0.
\end{align*}
Since $f_\Lambda$ is $1$-homogeneous, its derivative is $0$-homogeneous. Thus, if we replace $B$ by $\epsilon B$ and let $\epsilon\searrow$, we obtain
\begin{equation*}
d f_\Lambda(A)A\leq d f_\Lambda(B)A.
\end{equation*}
Moreover, the $1$-homogeneity of $f_\Lambda$ implies $d f_\Lambda(A)A=f_\Lambda(A)$, which settles the claim.
\end{proof}

\begin{theorem}[Duality formula]\label{thm:duality}
For $\rho_0,\rho_1\in \S_+(A)$ we have
\begin{align*}
\frac1  2\W_{\vec\Lambda,\epsilon}(\rho_0,\rho_1)^2&=\sup\{\tau(A(1)\rho_1)-\tau(A(0),\rho_0) : A\in \HJB_{\vec\Lambda,\epsilon}\}\\
&=\sup\{\tau(A(1)\rho_1)-\tau(A(0),\rho_0) : A\in \HJB_{\vec\Lambda,\epsilon}\cap C^\infty([0,1];\A)\}.
\end{align*}
\end{theorem}
\begin{proof}
The second inequality follows easily by mollifying. We will show the duality formula for Hamilton--Jacobi subsolutions in $H^1$. For this purpose we use the Rockefellar--Fenchel duality formula from Theorem \ref{thm:Rockefellar-Fenchel}.

Let $E$ be the real Banach space
\begin{equation*}
H^1([0,1];\Hh_\A)\times L^2([0,1];\Hh_{\A,\J}).
\end {equation*}
By the theory of linear ordinary differential equations, the map
\begin{equation*}
H^1([0,1];\Hh_\A)\to\Hh_\A\times L^2([0,1];\Hh_\A),\,A\mapsto (A(0),\dot A+\epsilon \L A)
\end{equation*}
is a linear isomorphism

Thus the dual space $E^\ast$ can be isomorphically identified with
\begin{align*}
\Hh_\A\times L^2([0,1];\Hh_\A)\times L^2([0,1];\Hh_{\A,\J})
\end{align*}
via the dual pairing
\begin{align*}
\langle (A,\bV),(B,C,\bW)\rangle&=\tau(A(0)B)+\int_0^1 \tau((\dot A(t)+\epsilon \L A(t))C(t))\,dt\\
&\quad+\int_0^1 \langle\bV(t),\bW(t)\rangle_{\H_{\A,\J}}\,dt.
\end{align*}
Define functionals $F,G\colon E\lra \IR\cup\{\infty\}$ by
\begin{align*}
F(A,\bV)&=\begin{cases}-\tau(A(1)\rho_1)+\tau(A(0)\rho_0)&\text{if }\bV=\nabla A,\\
\infty&\text{otherwise},
\end{cases}\\
G(A,\bV)&=\begin{cases}
0&\text{if }(A,\bV)\in \D,\\
\infty&\text{otherwise}.
\end{cases}
\end{align*}
Here
\begin{equation*}
\D=\{(A,\bV): \tau((\dot A(t)+\epsilon \L A(t))\rho)+\frac 1 2\norm{\bV(t)}_{\rho}^2\leq 0\text{ for all }t\in [0,1],\rho\in \S(\A)\}.
\end{equation*}
It is easy to see that $F$ and $G$ are convex. Moreover, for $A_0(t)=-t 1_\A$ and $\bV_0=0$ we have $\bV_0=\nabla A_0$, hence $F(A_0,\bV_0)=0$, and
\begin{equation*}
\tau((\dot A_0(t)+\epsilon \L A_0(t)) \rho)+\frac 1 2\norm{\bV_0(t)}_{\rho}^2=-1
\end{equation*}
for all $t\in[0,1],\;\rho\in\S(A)$, hence $G(A_0,\bV_0)=0$. Furthermore, $G$ is clearly continuous at $(A_0,\bV_0)$.

Moreover,
\begin{equation*}
\sup_{(A,\bV)\in E}(-F(A,\bV)-G(A,\bV))=\sup_{A\in\HJB_{\vec\Lambda,\epsilon}(\rho_0,\rho_1)}(\tau(A(1)\rho_1)-\tau(A(0)\rho_0)).
\end{equation*}

Let us calculate the Legendre transforms of $F$ and $G$, keeping in mind the identification of $E^\ast$. For $F$ we obtain
\begin{align*}
F^\ast(B,C,\bW)&=\sup_{(A,\bV)\in E}\bigg\{\langle (A,\bV),(B,C,\bW)\rangle-F(A,\bV)\bigg\}\\
&=\sup_{A}\bigg\{\tau(A(0)B)+\int_0^1 \tau((\dot A(t)+\epsilon \L A(t))\,C(t))\,dt\\
&\quad+\int_0^1 \langle\nabla A(t),\bW(t)\rangle_{\H_{A,\J}}\,dt+\tau(A(1)\rho_1)-\tau(A(0)\rho_0)\bigg\}.
\end{align*}
Since the last expression is homogeneous in $A$, we have $F^\ast(B,C,\bW)=\infty$ unless
\begin{align*
-\tau(A(1)\rho_1)+\tau(A(0)(\rho_0-B))&=\int_0^1 \tau((\dot A(t)+\epsilon \L A(t))\,C(t))\,dt\\
&\quad+\int_0^1 \langle \nabla A(t),\bW(t)\rangle_{\H_{\A,\J}}\,dt
\end{align*}
for $A\in H^1([0,1];\Hh_\A)$.

This implies $C_0=-(\rho_0-B)$ and $C_1=-\rho_1$ and
\begin{equation*}
\dot C(t)+\div \bW(t)=\epsilon \L^\dagger C(t).
\end{equation*}
Thus
\begin{equation*}
F^\ast(B,C,\bW)=\begin{cases}0&\text{if }(-C,-\bW)\in \CE_\epsilon^{\prime\prime}(\rho_0-B,\rho_1),\\
\infty&\text{otherwise}.
\end{cases}
\end{equation*}
Here $\CE_\epsilon^{\prime\prime}(\rho_0-B,\rho_1)$ denotes the set of all pairs $(X,\mathbf{U})\in H^1((0,1);\Hh_\A)\times L^2((0,1);\Hh_{\A,\J})$ satisfying $X(0)=\rho_0-B$, $X(1)=\rho_1$ and
\begin{equation*}
\dot X(t)+\div \bm{U}(t)=\epsilon \L^\dagger X(t).
\end{equation*}

The difference to the definitions of $\CE$ (or $\CE^\prime$) and $\CE^{\prime\prime}$ is that we do not make any positivity or normalization constraints. Note however that if $(X,\bm U)\in \CE^{\prime\prime}(\rho_0-B,\rho_1)$, then
\begin{equation*}
\frac{d}{dt}\tau(X(t))=\tau(\epsilon \L^\dagger X(t)-\div \bm U(t))=0
\end{equation*}
so that $\tau(X(t))=\tau(\rho_1)=1$ (and $\tau(B)=0$).

Now let us turn to the Legendre transform of $G$. We have
\begin{align*}
G^\ast(B,C,\bW)&=\sup_{(A,\bV)\in E}\bigg\lbrace \langle (A,\bV),(B,C,\bW)\rangle-G(A,\bV)\bigg\rbrace\\
&=\sup_{(A,\bV)\in \D}\bigg\lbrace \tau(A(0)B)+\int_0^1 \tau((\dot A(t)+\epsilon \L A(t))C(t))\,dt\\
&\qquad\qquad+\int_0^1\langle \bV(t),\bW(t)\rangle_{\H_{\A,\J}})\,dt\bigg\rbrace.
\end{align*}
Since $(A,\bV)\in \D$ implies $(A+X,\bV)\in \D$ for all $X\in \A$, we have $G^\ast(B,C,\bV)=\infty$ unless $B=0$. Moreover, it follows from the definition of $\D$ that $G^\ast(0,C,\bW)=\infty$ unless $C(t)\geq 0$ for a.e. $t\in [0,1]$.

For $B=0$ we have
\begin{align*}
G^\ast(0,C,\bW)&=\sup_{(A,\bV)\in\D}\left\lbrace\int_0^1 (\tau((\dot A(t)+\epsilon \L A(t))C(t))+\langle \bV(t),\bW(t)\rangle_{\H_{\A,\J}})\,dt\right\rbrace\\
&\leq \sup_{(A,\bV)\in\D}\bigg\lbrace -\int_0^1 \frac 1 2\norm{\bV(t)}_{C(t)}^2\,dt\\
&\qquad\qquad+\int_0^1\langle [C(t)]_{\vec\Lambda}^{1/2}\bV(t),[C(t)]_{\vec\Lambda}^{-1/2}\bW(t)\rangle_{\H_{\A,\J}}\,dt\bigg\rbrace\\
&\leq \frac 1 2\int_0^1\langle [C(t)]_{\vec\Lambda}^{-1} \bW(t),\bW(t)\rangle_{\H_{\A,\J}}\,dt.
\end{align*}
We will show next that the inequalities are in fact equalities. Let $C^\delta=C+\delta$
and $\bV^{\delta}(t)=[C(t)^{\delta}]^{-1}\bW(t)$. Moreover, let $f_j=f_{\Lambda_j}$ with the notation from Lemma \ref{lem:deriv_conv_fct}. Since
\begin{equation*}
\Hh_\A\to \IR,\,B\mapsto \sum_{j\in\J}\langle (df_{j}(C^{\delta}(t))B)\bV^{\delta}_j(t),\bV^{\delta}_j(t)\rangle_{\H_\A}
\end{equation*}
is a bounded linear map that depends continuously on $t$, there exists a unique continuous map $X^{\delta}\colon [0,1]\to \Hh_\A$ such that
\begin{equation*}
\tau(B X^{\delta}(t))=\sum_{j\in\J}\langle (df_{j}(C^{\delta}(t))B)\bV^{\delta}_j(t),\bV^{\delta}_j(t)\rangle_{\H_\A}
\end{equation*}
for every $B\in\Hh_\A$ and $t\in[0,1]$.

Let
\begin{equation*}
 A^{\delta}\colon [0,1]\to \A_h,\,A^{\delta}(t)=-\frac1  2\int_0^t X^{\delta}(s)\,ds.
\end{equation*}
We claim that $(A^{\delta},\bV^{\delta})\in \D$. Indeed,
\begin{align*}
\tau(\dot A^{\delta}(t)\rho)&=-\frac 1 2\sum_{j\in\J}\langle (df_{j}(C^{\delta}(t))\rho)\bV^{\delta}_j(t),\bV^{\delta}_j(t)\rangle_{\H_\A}\\
&\leq -\frac 1 2\sum_{j\in\J}\langle [\rho]_{\Lambda_j}\bV^{\delta}_j(t),\bV^{\delta}_j(t)\rangle_{\H_\A}\\
&=-\frac 1 2\norm{\bV^{\delta}(t)}_{\rho}^2,
\end{align*}
where the inequality follows from Lemma \ref{lem:deriv_conv_fct}. Note that we have equality for $\rho=C^{\delta}(t)$.

In particular, for $\rho=C(t)$ we obtain
\begin{equation*}
\tau(\dot A^{\delta}(t)C(t))+\langle \bV^{\delta}(t),\bW(t)\rangle_{\H_{\A,\J}}\leq \frac 1 2\langle [C(t)]_{\vec\Lambda}^{-1}\bW(t),\bW(t)\rangle_{\H_{\A,\J}}.
\end{equation*}
On the other hand,
\begin{align*}
\tau(\dot A^{\delta}(t)C(t))&=-\frac 1 2\sum_{j\in\J}\langle(d f_j(C^\delta(t))(C^\delta(t)-\delta))\bV_j^\delta(t),\bV_j^\delta(t)\rangle_{\H_A}\\
&\geq-\frac 1 2\langle [C^\delta(t)]_{\vec\Lambda}\bV^\delta(t),\bV^\delta(t)\rangle_{\H_{\A,\J}}+\frac 1 2\langle [\delta]_{\vec \Lambda}\bV^\delta(t),\bV^{\delta}(t)\rangle_{\H_{\A,\J}}\\
&\geq -\frac 1 2 \langle \bV^\delta(t),\bW(t)\rangle_{\H_{\A,\J}},
\end{align*}
where we again used Lemma \ref{lem:deriv_conv_fct} for the first inequality.

Put together, we have
\begin{align*}
\frac 1 2\langle [C^\delta(t)]_{\vec\Lambda}^{-1}\bW(t),\bW(t)\rangle_{\H_{\A,\J}}&\leq \tau(\dot A^{\delta}(t)C(t))+\langle \bV^{\delta}(t),\bW(t)\rangle_{\H_{\A,\J}}\\
&\leq \frac 1 2\langle [C(t)]_{\vec\Lambda}^{-1}\bW(t),\bW(t)\rangle_{\H_{\A,\J}},
\end{align*}
and
\begin{align*}
\int_0^1 (\tau(\dot A^\delta(t)C(t))+\langle \bV^\delta(t),\bW(t)\rangle_{\H_{\A,\J}})\,dt\to \frac 1 2\int_0^1 \langle [C(t)]_{\vec\Lambda}^{-1}\bW(t),\bW(t)\rangle_{\H_{\A,\J}}\,dt
\end{align*}
follows from the monotone convergence theorem.

Hence
\begin{equation*}
G^\ast(0,C,\bW)=\frac 1 2\int_0^1 \langle [C(t)]_{\vec\Lambda}^{-1}\bW(t),\bW(t)\rangle_{\H_{\A,\J}}\,dt
\end{equation*}
if $C(t)\geq 0$ for a.e. $t\in [0,1]$. Together with the formula for $F^\ast$, we obtain
\begin{align*}
&\quad\inf_{(B,C,\bW)\in E^\ast}(F^\ast(-B,-C,-\bW)+G^\ast(B,C,\bW))\\
&=\inf_{(\rho,\bV)\in \CE_\epsilon^\prime(\rho_0,\rho_1)}\frac 1 2 \int_0^1 \langle [\rho(t)]_{\vec\Lambda}^{-1}\bV(t),\bV(t)\rangle_{\H_{\A,\J}}\,dt\\
&=\frac 1 2 \W_{\vec\Lambda,\epsilon}^2(\rho_0,\rho_1),
\end{align*}
where the last equality follows from Proposition \ref{prop:relaxed_paths_CE}.

An application of the Rockefellar--Fenchel theorem yields the desired conclusion.
\end{proof}

\DeclareFieldFormat[article]{citetitle}{#1}
\DeclareFieldFormat[article]{title}{#1}

\printbibliograph

\end{document}


\section{Spherical Shape Analogies}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/analogies.jpg}
  \vspace{-5pt}
  \caption{We generate an output shape $B'$ that relates to the input $B$ in the same way as how the surface normal of a given primitive $A'$ relates to the surface normal of a sphere $A$. }
  \label{fig:analogies}
  \vspace{-10pt}
\end{figure}
Our main idea is to use surface normals to capture the style of 3D objects: if two shapes share a similar normal ``profile'', we consider them to exhibit the same geometric style.
Centered around this observation, as discussed in \refsec{intro}, we propose an analogy-based stylization method to translate the relationship between the normals of $A', A$ to create a stylized output shape $B'$ (see \reffig{analogies}).
Throughout the paper, we use green color to denote the target style shape $A'$, gray color to denote the input shape $B$, and blue color to denote the output stylized shape $B'$.

\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/analogies_gauss.jpg}
  \vspace{-5pt}
  \caption{Our algorithm defines the analogous relation based on the surface normals. We first map the normals of the style shape $N_{A'}$ to a unit sphere to obtain $\widetilde{N}_{A'}$ (top row), transfer the relationship between $N_A$ and $\widetilde{N}_{A'}$ to the input shape to obtain the target normals $T$ (middle row), then optimize the input shape $B$ so that the actual output normals are aligned with the target normal $T$ (bottom row). }
  \label{fig:analogies_gauss}
  \vspace{-10pt}
\end{figure}

Our algorithm consists of three simple steps, described in \reffig{analogies_gauss}: (1) we map the surface normal of $A'$ to a unit sphere $A$ in order to compute target normals on a sphere $\widetilde{N}_{A'}$, (2) we construct analogous target normals $T$ that relate to $N_{B}$ the same way $\widetilde{N}_{A'}$ relate to $N_{A}$, (3) we take $B, T$ as inputs and generate the stylized shape $B'$ whose normals approximate $T$ via optimization.

\subsection{Generating $\widetilde{N}_{A'}$}\label{sec:generateNS}
Depending on the provided style shape $A'$ or user preferences, we consider three ways to get a set of target normals on a sphere $\widetilde{N}_{A'}$.

\textit{1. Closest normals.}\
The simplest case is when the style shape $A'$ is a simple convex shape with only few distinct face normals (e.g., icosahedron). We compute $\widetilde{N}_{A'}$ simply via snapping the normals of the sphere $N_{A}$ to the nearest normal in the style shape $N_{A'}$.

\textit{2. Spherical parameterization.}\
For a generic genus-0 shape (e.g., smooth or concave), we compute its parameterization to a sphere using, for example, conformalized mean curvature flow \cite{KazhdanSB12}. Then $\widetilde{N}_{A'}$ can be computed from the spherical parameterization.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/meanCurvatureFlow.jpg}
  \vspace{-5pt}
  \caption{Given a style shape $A'$, we run mean curvature flow \cite{KazhdanSB12} to map the normals of style shape $N_{A'}$ to a sphere as $\widetilde{N}_{A'}$.}
  \label{fig:meanCurvatureFlow}
  \vspace{-10pt}
\end{figure}

\textit{3. Normal Capture.}\
If one desires more control, one can manually specify $\widetilde{N}_{A'}$, in the spirit of how \emph{MatCap} (material capture \cite{SloanMGG01}) is used in rendering. We can then skip the first step in \reffig{analogies_gauss} and move on to the second step using the user-provided $\widetilde{N}_{A'}$.

\subsection{Generating $T$}\label{sec:generateT}
\begin{wrapfigure}[6]{r}{1.33in}
	\raggedleft
  \vspace{-12pt}

	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.48in, trim={6mm 0mm -1mm 0mm}]{figures/gaussmap.jpg}
	\label{fig:gaussmap}
\end{wrapfigure}
Generating target normals $T$ on the input shape $B$ using analogy requires the correspondences between $A$, $B$. We compute the map using the \emph{Gauss map}, leveraging the fact that our $A$
is always a unit sphere (see the inset, where we use colors to visualize the correspondences).
\update{Specifically, the unit normal vector of each element (e.g., vertex or face) on the input shape $B$ can be equivalently interpreted as a point on the unit sphere $A$. Thus, we can easily map signals from $A$ back to $B$.}
Once the correspondences are obtained via the normals of input shape $N_{B}$, we can trivially compute $T$ by ``pasting'' $\widetilde{N}_{A'}$ on top of $B$.

\subsection{Generating $B'$}
After obtaining a set of target normals $T = \{ {\vecFont{t}}_k \}$ for each vertex $k$, our goal is to obtain a deformed output shape $B'$ whose surface normals approximate $T$.
Let ${\matFont{V}}$ be a matrix of vertex locations with size $|{\matFont{V}}|$-by-3 and ${\matFont{F}}$ be the face list with size $|{\matFont{F}}|$-by-3 of the input shape $B$. Our output shape $B'$ is a deformed version of the input shape and we use $\widetilde{{\matFont{V}}}$ to denote the $|{\matFont{V}}|$-by-3 matrix of the deformed vertex locations. We formulate the normal-driven deformation as an energy optimization in the following form:
\begin{align}
  \min_{\widetilde{{\matFont{V}}}} \sum_{k \in {\matFont{V}}} E_R({\vecFont{v}}_k,{\vecFont{u}}_k) + \lambda a_k \| \hat{n}_k(\widetilde{{\matFont{V}}}) - {\vecFont{t}}_k \|_2^2,
\end{align}
where $E_R$ denotes a regularization energy to preserve the details of the input mesh, and the second part measures the squared distance from the output unit surface normal $\hat{{\vecFont{n}}}_k(\widetilde{{\matFont{V}}})$ to the target output normal ${\vecFont{t}}_k$ at vertex $k$. We use $a_k$ to denote the Voronoi area of the vertex $k$, $\lambda$ is a weighting parameter to control the balance between the two terms, and ${\vecFont{v}}_k$ (${\vecFont{u}}_k$) is the input (output) location of vertex $k$.
In \reffig{lambda}, we can observe that using a small $\lambda$, the method preserves the input shape $B$. Using a bigger $\lambda$, the method favors in deforming the shape more into the style of $A'$.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/lambda_more.jpg}
  \vspace{-5pt}
  \caption{\update{The $\lambda$ parameter in \refequ{arapNormalDriven} controls the balance between preserving the original shape and satisfying the desired style. We show different stylization results with increasing $\lambda$. \copyright Spiral Light Bulb (top) by benglish under CC BY-SA.}}
  \label{fig:lambda}
  \vspace{-10pt}
\end{figure}
The choice of $E_R$ depends on the user's intent. One can apply different regularizations to obtain different results.
For the purposes of this exposition, we introduce our optimization based on {\textsc{arap}}\xspace regularization in \refsec{normalDriven_arap}. We discuss how to extend to other regularizations in \refsec{regularizations}.

\subsection{Normal-Driven Optimization with {\textsc{arap}}\xspace}\label{sec:normalDriven_arap}
We use ${\vecFont{e}}_{ij} \coloneqq {\vecFont{v}}_j - {\vecFont{v}}_i \in \mathbb{R}^3$ to denote the edge vector between vertices $i, j$ on the original mesh, and ${{\vecFont{e}}'}_{ij} \coloneqq \u_j - \u_i$ for the edge vectors on the deformed mesh. We can write down the energy that uses {\textsc{arap}}\xspace regularization as
\begin{align}\label{equ:arapOutputNormal}
  \min_{\widetilde{{\matFont{V}}},{\matFont{R}}}\ \sum_{k \in  V} \underbrace{\sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2}_{E_{\textsc{arap}}\xspace} + \lambda a_k \| \hat{n}_k(\widetilde{{\matFont{V}}}) - {\vecFont{t}}_k \|_2^2,
\end{align}

\begin{wrapfigure}[6]{r}{0.9in}
	\raggedleft

	\vspace{-14pt}
	\hspace*{-0.9\columnsep}
	\includegraphics[width=1.25in, trim={5mm 0mm 1mm 0mm}]{figures/spokes_rims.jpg}
\end{wrapfigure}
We use $\mathcal{N}_k$ to denote the edge vectors of the \emph{spokes and rims} at vertex $k$ (see the inset) \cite{ChaoPSS10}, ${\matFont{R}}_k \in \text{SO(3)}$ to denote a 3-by-3 rotation matrix defined on $k$, and $w_{ij}$ is the cotangent weight of edge $i,j$ \cite{PinkallP93}.
However, this energy is difficult to optimize because the term $\hat{n}_k(\widetilde{{\matFont{V}}})$ is non-linear in $\widetilde{{\matFont{V}}}$.

We adapt the observation made in \cite{LiuCubic19} that the space of unit vectors can be captured by rotations. Thus, we can perform a change of variables by replacing $\hat{n}_k(\widetilde{{\matFont{V}}})$ with the \emph{rotated} unit normal of the input mesh ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ as
\begin{align} \label{equ:arapNormalDriven}
  \boxed{
  \min_{\widetilde{{\matFont{V}}},{\matFont{R}}}\ \sum_{k \in V} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2  + \lambda a_k \| {\matFont{R}}_k \hat{{\vecFont{n}}}_k  - {\vecFont{t}}_k \|_2^2,}
\end{align}
where $\hat{{\vecFont{n}}}_k$ is the $k$th unit vertex normal of the input mesh computed via area-weighted average of face normals, which is constant throughout the optimization. This ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ can be perceived as an approximation of the area-weighted vertex normals of the output mesh $\hat{n}_k(\widetilde{{\matFont{V}}})$. In \reffig{RN_difference}, we visualize the difference between the output normals $\hat{n}_k(\widetilde{{\matFont{V}}})$ and the rotated input normals ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$. We can notice that ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ is a decent approximation of the output vertex normals computed via area-weighted average.
\update{We can observe that error tend to concentrate on high-curvature regions because discrete vertex normals are less accurate along on those regions and the {\textsc{arap}}\xspace regularization encourages smooth deformation.}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/RN_difference.jpg}
  \vspace{-5pt}
  \caption{Empirically, we show that rotated input normals ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ is a good approximation of the area-weighted output vertex normals $\hat{n}_k(\widetilde{{\matFont{V}}})$. We can observe that the error mostly occurs on the high-curvature regions (right). \copyright Proto Paradigm under CC BY.}
  \label{fig:RN_difference}
  \vspace{-10pt}
\end{figure}
This change of variables allows us to solve for ${\matFont{R}}_k$s in parallel and make this energy quadratic in $\widetilde{{\matFont{V}}}$. In addition, the fact that ${\matFont{R}}_k$ is shared across the {\textsc{arap}}\xspace term and the normal term enables us to jointly consider both the regularization and the normal terms when obtaining the deformed vertex locations $\widetilde{{\matFont{V}}}$.

We minimize this energy via the local/global strategy \cite{SorkineA07}, where the local step involves solving a set of small \emph{Orthogonal Procrustes} problems and the global step amounts to a linear solve. For the sake of reproducibility, we reiterate the local-global steps for our energy in \refapp{localStep}, \refappnum{globalStep}.
Non-linear methods, such as Newton’s method, could be applied to our scenario. It is however far slower than the local-global optimization since a single iteration of the Newton's method could be more expensive than 100 iterations of the local-global iterations (see \cite{LiuBK17}). Thus, it is less suitable for our interactive applications.
Further accelerating our solver using other optimization methods (e.g., \cite{KovalskyGL16, PengDZGQL18, ZhuBK18}) should be possible, but is left as future work.

\section{Introduction}\label{sec:intro}

Analogies of the form $A:A'\dblcolon B:B'$ is a reasoning process that conveys \emph{$A$ is to $A'$ as $B$ is to $B'$}. This formulation has become a core technique for creating artistic 2D digital content, such as image analogies \cite{HertzmannJOCS01} in Photoshop \cite{photoshop} for image stylization and the Lit Sphere \cite{SloanMGG01} (a.k.a. MatCap) in ZBrush \cite{zbrush} for non-photorealistic renderings.
However, leveraging analogies to stylize 3D geometry is still at a preliminary stage because defining the analogy relationship on surface meshes requires dealing with irregular discretizations, curved metrics, and different topologies.

In this paper, we introduce a step towards a more general 3D shape analogies, named \emph{spherical shape analogies}. We consider a specific case where $A$ is a unit sphere. This restriction enables us to operate on an input mesh $B$ with arbitrary topologies, boundaries, and geometric complexity.
While not fully general, because $A$ is restricted to be a sphere, we demonstrate that this formulation can immediately achieve different geometric styles within a single framework. In \reffig{teaser}, we show that by providing different target style shapes $A'$ to the algorithm, we can turn the input shape $B$ into different styles.
In addition to stylization, our method can encompass many existing applications, such as developable surface approximation and PolyCube deformation.

One key observation in our spherical shape analogies is that the surface normal is an effective instrument to capture geometric styles. Thus, we define the analogy relationship based on normals: we optimize a stylized shape $B'$ such that the relationship between the surface normals of $B$ and $B'$ is the same as the relationship between the surface normals of $A$ and $A'$

We realize this by casting it as a simple and effective normal-driven shape optimization problem which aims at deforming the input shape towards a set of desired normals. However, such an optimization problem is often difficult due to the nonlinearity of unit normals. We draw inspiration from previous works and apply a change of variables to accelerate the computation: instead of directly optimizing the vertex positions, we optimize a set of rotations that rotate the normals of the input mesh to the set of desired normals.
Our simple formulation with the change of variables results in a generic stylization algorithm that runs at interactive rates.

\section{Extensions \& Analysis}\label{sec:analysis}
In this section, we introduce its extensions to different regularizations and how to handle cases where target normals $T(B')$ are a function of output geometry.

\subsection{Different Regularizations}\label{sec:regularizations}
In addition to $E_{\textsc{arap}}\xspace$, the normal-driven optimization supports different regularization energies for different modeling intents.
One could use {\textsc{arap}}\xspace when the goal is to produce a smooth deformation that preserves surface details.
If one wants to produce a non-smooth deformation (e.g., sharp creases) while preserving local rigidity, one could instead use a \emph{face-only} {\textsc{arap}}\xspace energy $E_{\textsc{farap}}\xspace$ \update{discussed in \cite{ZhaoG16a, LeviG15}} which consists of only the membrane term.
If one is interested in preserving the textures and allowing local scaling, one could use an \emph{as-conformal-as-possible} energy $E_{\textsc{acap}}\xspace$ \cite{BouazizDSWP12}.

\begin{wrapfigure}[4]{r}{0.53in}
	\raggedleft
	\vspace{-16pt}

	\hspace*{-0.9\columnsep}
	\includegraphics[width=0.85in, trim={3mm 0mm 1mm 0mm}]{figures/face_only.jpg}
	\label{fig:stretchElement}
\end{wrapfigure}
\textit{Face-only {\textsc{arap}}\xspace.}\
The core idea is to remove the bending term from {\textsc{arap}}\xspace and only measure the membrane term \cite{TerzopoulosPBF87}, so that two adjacent triangles can bend freely. We achieve this by applying the idea from \cite{ZhaoG16a, LeviG15} which only measures the {\textsc{arap}}\xspace energy over the three edge vectors of a face $f_k$ (see the inset), instead of the spokes and rims $\mathcal{N}_k$. Precisely, we can write this ``face-only'' {\textsc{arap}}\xspace regularization $E_{\textsc{farap}}\xspace$ as
\begin{align}\label{equ:farap}
  E_{\textsc{farap}}\xspace(\widetilde{{\matFont{V}}}, {\matFont{R}}) =  \sum_{k \in F} \sum_{i,j \in f_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2,
\end{align}

\textit{As-conformal-as-possible.}\
If the goal is to create novel geometric details, it is crucial to allow non-rigid deformations. However, an arbitrary deformation may lead to undesirable behaviors, such as badly shaped triangles. Thus constraining the angle preservation, a.k.a. conformality, will be a suitable regularization.
Specifically, we use the {\textsc{acap}}\xspace energy $E_{\textsc{acap}}\xspace$ in \cite{BouazizDSWP12} as our regularization
\begin{align}\label{equ:acap}
  E_{\textsc{acap}}\xspace(\widetilde{{\matFont{V}}}, {\matFont{R}}, {\vecFont{s}}) =  \sum_{k \in F} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| s_k {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2,
\end{align}
where $s_k$ is a scalar representing the scaling of local patch. One can compute the optimal $s_k$ analytically via the method by Sch{\"o}nemann et al. \cite{schonemann1970fitting} (see \refapp{regularizations}).

Deploying these regularizations $E_{\textsc{farap}}\xspace, E_{\textsc{acap}}\xspace$ requires only a few changes in the optimization steps. Deploying $E_{\textsc{farap}}\xspace$ only involves changing the incidence matrix. Deploying $E_{\textsc{acap}}\xspace$ only requires adding one more line of code in the local step to solve an \emph{isotropic orthogonal Procrustes} problem \cite{schonemann1970fitting}. We detail such changes in \refapp{regularizations}. In \reffig{regularizations}, we apply the same deformation to a sheet but with different regularizations. We can perceive that different regularizations favor drastically different solutions.

Our framework allows one to easily plug-and-play different regularizations. Specifically, we use $E_{\textsc{arap}}\xspace$ for applications that favor smooth deformation (e.g., \reffig{primStyles}), $E_{\textsc{farap}}\xspace$ for creating sharp creases (\reffig{generalized_PolyCube}, \reffignum{developable}), and $E_{\textsc{acap}}\xspace$ when one wants to manipulate geometric details such as in \reffig{synthesis}.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/regularizations.jpg}
  \vspace{-5pt}
  \caption{Different regularizations favor different behaviors. Given a sheet (gray), we pull up the center part (central blue dots) and shrink the boundary (blue dots on the boundary), then we minimize each regularization energy to determine the unconstrained vertices. We can observe that $E_{\textsc{arap}}\xspace$ favors rigid and smooth interpolation, $E_{\textsc{farap}}\xspace$ favors sharp bending between triangles, and $E_{\textsc{acap}}\xspace$ favors to preserve angles while allowing local scaling. }
  \label{fig:regularizations}
  \vspace{-10pt}
\end{figure}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/update_t.jpg}
  \vspace{-5pt}
  \caption{Setting the target normal $T$ as a constant or treating it as a function of the output mesh $T(B')$ leads to different local minima. In many cases (left pair), both options lead to similar looking results. But setting $T$ as a constant may result in an undesirable local minimum in some cases (right), such as the ears of the bunny.}
  \label{fig:update_t}
  \vspace{-10pt}
\end{figure}

\subsection{Dynamic Target Normals}
Our method converges to a local minimum.
Empirically, we observe that treating the target normal $T$ as a constant throughout the optimization may work fine perceptually in many cases (see the left pair in \reffig{update_t}). However, constant $T$ may lead to an undesirable local minimum due to a sub-optimal assignment of $T$ (see the right pair in \reffig{update_t}).
Inspired by Projective Dynamics \cite{BouazizMLKP14}, a simple solution to avoid such local minima is to treat $T$ as a function of $B'$ ($N_{B'}$ specifically), and update $T$ at every iteration.
We summarize the pseudo code in \refalg{normalDriven}. If $T$ is a constant throughout the optimization, one can simply skip the optional step at line 9.
\IncMargin{1em}
\begin{algorithm}[t]
    \SetKwInOut{Input}{Input}
    \SetKwInOut{Output}{Output}
    \caption{Normal-driven optimization}
    \label{alg:normalDriven}
    \Indentp{-1em}
        \Input{$\,$A triangle mesh ${\matFont{V}}, {\matFont{F}}$ and a weight $\lambda$}
        \Output{$\,$Deformed vertex positions $\widetilde{{\matFont{V}}}$}
        \BlankLine
    \Indentp{1em}
        compute $\widetilde{N}_{A'}$ \hfill \text{\color{iglGreen} // step 1, Sec. 3.1}\\
        compute $T$ from $\hat{n}({\matFont{V}})$ \hfill \text{\color{iglGreen} // step 2, Sec. 3.2}\\
        $\hat{{\vecFont{n}}} \leftarrow \hat{n}({\matFont{V}})$ \hfill \text{\color{iglGreen} // compute input surface normals}\\
        ${\matFont{Q}}, {\matFont{K}} \leftarrow precompute({\matFont{V}},{\matFont{F}})$ \hfill \text{\color{iglGreen} // see App. B}\\
        $\widetilde{{\matFont{V}}} \leftarrow {\matFont{V}}$\\
        \While{\textit{not converge}}{
          ${\matFont{R}} \leftarrow local\text{-}step\,(\widetilde{{\matFont{V}}}, \hat{{\vecFont{n}}}, T, \lambda)$  \hfill \text{\color{iglGreen} // App. A}\\
          $\widetilde{{\matFont{V}}} \leftarrow global\text{-}step\,({\matFont{R}}, {\matFont{Q}},{\matFont{K}})$ \hfill \text{\color{iglGreen} // App. B}\\
          compute $T$ from $\hat{n}(\widetilde{{\matFont{V}}})$ \hfill \text{\color{iglGreen} // (optional) for dynamic $T$}\\
        }
\end{algorithm}
\DecMargin{1em}

\begin{wrapfigure}[7]{r}{1.33in}
	\raggedleft

  \vspace{-12pt}

	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.55in, trim={6mm 0mm -1mm 0mm}]{figures/convergence.jpg}
	\label{fig:convergence}
\end{wrapfigure}
In terms of convergence, in the case where $T$ is constant, the convergence behaves the same as the original {\textsc{arap}}\xspace \cite{SorkineA07}, where the energy decreases monotonically. In the case where $T$ is dependent to $B'$, we do not guarantee a monotonic decrease in energy, but the optimization still converges in our experiments. In the inset, we visualize the convergence plot for examples in \reffig{update_t}.

\begin{wrapfigure}[7]{r}{1.33in}
	\raggedleft
  \vspace{-12pt}

	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.55in, trim={6mm 0mm -1mm 0mm}]{figures/runtime.jpg}
	\label{fig:runtime}
\end{wrapfigure}
We implement our algorithm in C++ with Eigen \cite{eigenweb} and evaluate our method on a MacBook Pro with an Intel i5 2.3GHz processor.
Our method runs 24 iterations per second for a mesh with around 20k vertices. We report a complete picture of our runtime in the inset. The local step will be the computation bottleneck for meshes with less than 20k vertices, but further acceleration can be achieved via the method by Zhang et al. \cite{ZhangFast2021}.
Typically, within the first 10 iterations, our method can achieve a visually similar result compared to the converged solution. This property enables us to build an interactive tool for users to play with different style shapes $A'$ or artistic controls.

\section{Applications}
The major benefit of our analogy-based stylization method is that one can plug-and-play different style shapes to obtain different results. \update{When one provides} convex primitives with few distinct face normals, we can simply use the method discussed in \refsec{generateNS} to turn an input shape into the style of the primitive (see \reffig{primStyles}, \reffignum{primStyles_bd}). In \reffig{quan_normal}, we also quantitatively show that our method can effectively reduce the difference between mesh normals and the normals of a primitive.
If the provided style shape is smooth or non-convex, where the simple closest normal may fail to capture the style, one could use a spherical parameterization described in \refsec{generateNS} to achieve the stylization.
If desiring more user controls, one could ``paint'' the desired surface normals on a unit sphere (see \refsec{generateNS}), and then transfer the style of the painted normals directly to the input (see \reffig{manual_gauss})

\subsection{PolyCube Deformation}
If one is interested in PolyCube maps \cite{TariniHCM04}, we can adapt normal driven editing to create PolyCube maps,following the observation in \cite{ZhaoLLZXG17}. Specifically, we need to use a cube as a style shape and move the pre-computation step in \refalg{normalDriven} to the optimization loop.
Moving the pre-computation in the loop would no longer preserve the original details, which is desirable for creating PolyCube shapes. This modification may also lead to badly shaped triangle.
When these faces appear, a quick solution is to move the vertex towards the 1-ring average by a small amount to improve triangle quality.
For the sake of comparison, we use the same PolyCube segmentation as in \cite{FuBL16} and show that we can achieve comparable results in \reffig{PolyCube_comparison}.
We can further generalize the PolyCube map to other polygonal boxes by specifying non-cube normals (see \reffig{generalized_PolyCube}).

\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/generalized_PolyCube.jpg}
  \vspace{-5pt}
  \caption{\update{By using} different sets of normals, we can generalize the PolyCube method (left) to create polygonal boxed maps.}
  \label{fig:generalized_PolyCube}

  \vspace{-10pt}
\end{figure}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/control_developable.jpg}
  \vspace{-5pt}
  \caption{Stein et al. \cite{SteinGC18} control the patches on the developable surfaces via remeshing the input. We, instead, can control the amount of creases (middle, right) by tuning a single parameter (see \refapp{developable}).}

  \label{fig:control_developable}

  \vspace{-10pt}
\end{figure}

\subsection{Developable Surface Approximation}
So far we have only considered an explicit shape or a set of painted normals as our style shape. Here we further extend our method to support an energy that describes a certain style. In particular, we consider the target normal $T$ is computed via an optimization
\begin{align}
  T = \argmin_T f(B'),
\end{align}
and, similar to the case where $T$ is dependent to $B'$, we update $T$ at every iteration in the local/global solve.

\begin{wrapfigure}[7]{r}{1.33in}
	\raggedleft
    \vspace{-14pt}
	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.48in, trim={6mm 0mm -1mm 0mm}]{figures/K_hist}
	\label{fig:K_hist}
\end{wrapfigure}
We evaluate this extension via setting $f$ to be the discrete developability energy proposed in \cite{SteinGC18}, with details provided in \refapp{developable}. Compared to the original method, our approach contains a regularization term in addition to the developable energy, thus our optimization requires no remeshing and results in the faster optimization (see \reffig{developable_comaprison}).
In \reffig{control_developable}, we further show that our framework enables one to control the number of creases in the results.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/developable.jpg}
  \vspace{-5pt}
  \caption{Our normal driven editing can be used to create many piece-wise developable surfaces (blue). Our method requires no remeshing and is fast enough for interactive modeling. \copyright cerberus333 (third) under CC BY-NC.}
  \label{fig:developable}
  \vspace{-10pt}
\end{figure}
With our framework one can interactively create a variety of piece-wise developable shapes (see \reffig{developable}).
In \reffig{K_developable}, we evaluate our results by visualizing the discrete Gaussian curvature before and after running our developable flow. We can observe that the Gaussian curvature concentrates along the creases and results in a piece-wise developable surface. In the inset, we quantitatively demonstrate that our method effectively increases the developability of the mesh in \reffig{K_developable}.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/K_developable.jpg}
  \vspace{-5pt}
  \caption{We use our normal driven editing to deform the input shape (gray) into a piece-wise developable approximation (blue). In the bottom row, we visualize the Gaussian curvature concentrates on the creases after the deformation, leading to a piece-wise developable shape. \copyright Oliver Laric under CC BY-NC-SA.}
  \label{fig:K_developable}
  \vspace{-10pt}
\end{figure}

\section{Related Work}
Our work shares similar motivations to computer-assisted image stylization pioneered by Haeberli \cite{Haeberli90}. But since our outputs are stylized 3D geometries, we focus the discussion on geometric stylization and geometric deformation methods.

\subsubsection*{Analogy-based Geometric Stylization}
Many generative models have been proposed for creating stylized 3D objects, such as collage art \cite{gal20073d, Theobalt2007Collage}, manga style \cite{ShenLHC12}, cubic style \cite{LiuCubic19}, and neuronal homunculus \cite{ReinertRS12}. However, these methods are tailor-made for only a specific style.

\emph{Analogy} $A:A'\dblcolon B:B'$ is a powerful idea to achieve different stylization results within a single framework.
This idea has inspired several design tools for images \cite{HertzmannJOCS01}, non-photorealistic renderings \cite{SloanMGG01,Shechtman16}, and curves \cite{HertzmannOCS02}.
\update{Beyond 2D data, the idea of analogy has also been used for transferring 3D geometric details from one shape to another.
We omit the discussion on methods that are not based on analogies, such as mesh cloning \cite{ZhouHWTDGS06, TakayamaSSIBS11} and geometric learning \cite{LiuKCAJ20, HertzHGC20, WangAKCS20, chen2021decor, li_cvpr21}, and focus on analogy-based techniques.}
Ma et al. \cite{MaHSKW14} propose a method for 3D style transfer based on patch-based assembly. However, their method cannot handle free-form deformations and requires the source and the exemplar shape to share a similar structure in order to compute high-quality correspondences.
\update{Bhat et al. \cite{BhatIT04} propose a voxel-based texture synthesis method for transferring geometric details encoded in the volumetric grid. Berkiten et al. \cite{BerkitenHSMLR17} use metric learning for details represented as displacement maps. These methods are designed for high-frequency details (e.g., wrinkles on the surface).} In contrast, our spherical shape analogies focuses on larger scale free-form deformations.
Albeit limited --- in our analogies $A$ is restricted to the unit sphere --- our method enables a first step in this exciting direction.

\subsubsection*{Surface Normals in Shape Deformation}
A key insight of our spherical shape analogies is to leverage surface normals to capture geometric styles.
The surface normal is a fundamental geometric quantity and is ubiquitous in geometry processing.
A representative example is in the \emph{PolyCube} deformation \cite{TariniHCM04} where the goal is to optimize surface normals to be axis-aligned.
Gregson et al. \cite{GregsonSZ11} and Zhao et al. \cite{ZhaoLLZXG17} use the closest rotation from the surface normal to an axis-aligned direction to drive the PolyCube deformation. Huang et al. \cite{HuangJSTBD14} and Fu et al. \cite{FuBL16} propose to minimize energies defined on normals to create PolyCube shapes.
In architectural geometry design, surface normals are a main ingredient to characterize polygon meshes with planar faces. The methods proposed by Deng et al. \cite{DengPW11} and Poranne et al. \cite{PoranneOG13} utilize normals to formulate a \emph{distance-from-plane} constraint to encourage planarity. Tang et al. \cite{TangSGWP14} use the dot product between a face normal and its adjacent edge vectors to determine whether the vertices of a polygon are coplanar.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/silvia.jpg}
  \vspace{-5pt}
  \caption{Sell{\'{a}}n et al. \cite{Sellan2020} propose a technique to make 2D heightfields developable (purple). In contrast, our method can create developable approximations for surface meshes in 3D (blue). }
  \label{fig:silvia}
  \vspace{-10pt}
\end{figure}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/developable_comaprison.jpg}
  \vspace{-5pt}
  \caption{Compared to the method proposed by Stein et al. \cite{SteinGC18} for creating developable approximations (left), our method can create visually comparable results (right) with significant speed-ups.}
  \label{fig:developable_comaprison}
  \vspace{-10pt}
\end{figure}
Characterizing whether a mesh can be flattened to 2D without stretching or shearing, a.k.a. \emph{developability}, also relies on surface normals. Stein et al. \cite{SteinGC18} characterize discrete developability based on the 1-ring face normals, and propose an algorithm to compute piecewise developable surfaces. Sell{\'{a}}n et al. \cite{Sellan2020} reformulate the developable energy into a convex semidefinite program for finding piecewise developable heightfields.
In addition to these examples, deforming shapes into the cubic style \cite{LiuCubic19, Fumero0R20}, constructing shape abstractions \cite{alexa2021polycover}, surface parameterization \cite{ZhaoSLZYLWGG20}, \update{and interactive mesh editing \cite{YuZXSBGS04, SorkineCLARS04}} are all related to surface normals. Many more examples can be found in the design of geometric filters, such as the Guided filter \cite{ZhangDZBL15}, the Shock filter \cite{PradaK15}, the Bilateral normal filter \cite{ZhengFAT11}, and the Total Variation mesh denoising \cite{ZhangWZD15}.

Our method can be adapted to these normal-based deformations. Compared to the PolyCube method \cite{FuBL16}, we achieve comparable quality (see \reffig{PolyCube_comparison}), but we can further generalize to polytopes (see \reffig{generalized_PolyCube}).
Compared to \cite{LiuCubic19} in cubic stylization (see \reffig{cubic_comparison}), we can achieve similar performance, but we can further generalize to many styles other than the cubic style (see \reffig{teaser}).
In developable surface approximation, in contrast to the method by Sell{\'{a}}n et al. \cite{Sellan2020}, our method can be applied to surface triangle meshes (see \reffig{silvia}) and is significantly faster than the method by Stein et al. \cite{SteinGC18} (see \reffig{developable_comaprison}).

\subsubsection*{Shape Deformation}
Our geometric stylization method can also be perceived as \update{a type of shape deformation method}. We share technical similarities with methods that deform a shape while addressing given modeling constraints.
A common choice is to minimize the \emph{as-rigid-as-possible} ({\textsc{arap}}\xspace) energy \cite{SorkineA07,IgarashiMH05, ChaoPSS10} while satisfying the constraints. This {\textsc{arap}}\xspace energy measures the rigidity of local surface patches and favors detail-preserving smooth deformations.
In the case where locally rigid deformations are too constrained, the conformal energy \cite{CranePS11, VaxmanMW15} which preserves angles is commonly used.  In contrast to {\textsc{arap}}\xspace, the conformal energy often triggers larger deformations as it allows both local uniform scaling and rigid transformations.
In addition to mesh deformations, similar energies have also been used for parameterization \cite{LiuZXGG08}, shape optimization \cite{BouazizDSWP12}, and simulating mass-spring systems \cite{LiuBOK13}.
The {\textsc{arap}}\xspace and conformal energies are also commonly used as regularization terms in mesh optimization problems, such as reconstruction \cite{IzadiF14}, surface registration \cite{HuangAWG08, YoshiyasuMYK14}, PolyCube construction \cite{HuangJSTBD14}, and surface stylization \cite{LiuCubic19}.
Their popularity comes from the property that they favor smooth deformations and are amenable to fast optimizations.
For the same reasons, we also use these as our regularization energies for interactive modeling tasks (see \reffig{regularizations}).

\section{Limitations \& Future Work}\label{sec:limitations}
Our method draws inspiration from Projective Dynamics \cite{BouazizMLKP14} to handle the case where target normals $T$ are a function of output shape $B'$ (e.g., \reffig{update_t}, \reffignum{developable}). Although being fast and suitable for our intended interactive applications, it often struggles to converge to a highly accurate solution.
Extending our optimization to, for example, Newton's method would be desirable for applications that desire highly accurate solutions.

Our approach is restricted to a sphere as our reference shape $A$, and uses the Gauss map to determine the correspondences between $A$ and the input $B$. As the Gauss map purely relies on surface normals to determine the map, the resulting map is ignorant to area distortion.
This characteristic is beneficial to handle input shapes $B$ that are very different (e.g., different genus) from a sphere because in these cases it is challenging to obtain a map with low area distortion.
However, the price we have to pay is that we cannot support structured
\begin{wrapfigure}[8]{r}{0.62in}
	\raggedleft

	\vspace{-6pt}
	\hspace*{-0.7\columnsep}
	\includegraphics[width=0.75in, trim={4mm 0mm 1mm 0mm}]{figures/printedHead_cut}
\end{wrapfigure}
and high-frequency patterns (e.g., geometric texture synthesis). Thus, if one is interested in stylizing shapes with detailed textures, we suggest to first synthesize target normals on the surface directly \cite{WieLKT09} then perform the normal-driven optimization (\refsec{normalDriven_arap}).
In \reffig{synthesis} we demonstrate this alternative by unbaking an existing normal map for manufacturing purposes (see the inset) and synthesizing normal textures from an image.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/synthesis.jpg}
  \vspace{-5pt}
  \caption{If one is interested in creating high-frequency geometric textures, we recommend to compute target normals via texture synthesis and then optimize the geometry via the normal-driven optimization. We demonstrate an example of ``unbaking'' normal maps (left) and an example of geometric texture synthesis (right).}

  \label{fig:synthesis}
  \vspace{-10pt}
\end{figure}

Our method currently supports manifold triangle meshes.
Extending to non-manifold meshes, polygon meshes, volumetric meshes, and point clouds could be beneficial to handle real-world geometric data.
Not every shape or normal capture sphere is valid to serve as the style shape of our algorithm. Discovering the validity of a style shape is important to understand the behavior of these novel modeling methods.
Removing the assumption about the source shape being a sphere could lead to a more general analogy-based shape editing.
Based on the observation that surface normals are a promising geometric quantity to capture the style of a shape. Developing a better categorization of styles based on normals or exploring learning-based techniques on normals (instead of vertices) could lead to novel stylization methods.

\section*{Acknowledgements}
  Our research is funded in part by NSERC Discovery (RGPIN2017–05235, RGPAS–2017–507938), New Frontiers of Research Fund (NFRFE–201), the Ontario Early Research Award program, the Canada Research Chairs Program, the Fields Centre for Quantitative Analysis and Modelling and gifts by Adobe Systems, Autodesk and MESH Inc.
  %
  We thank Sheldon Andrews, Abhishek Madan, Silvia Sell{\'{a}}n, Oded Stein, Li-Yi Wei for helps on experiments. We thank members of Dynamic Graphics Project at the University of Toronto; Sarah Kushner, Abhishek Madan, Silvia Sell{\'{a}}n, Let{\'{i}}cia Mattos Da Silva, Towaki Takikawa for proofreading; John Hancock for the IT support. We thank all the artists for sharing a rich variety of 3D models.
  %
\section{Local Step with $E_{\textsc{arap}}\xspace$} \label{app:localStep}
Given a fixed $\widetilde{{\matFont{V}}}$, we obtain the optimal rotation for each vertex $k$ by solving the following minimization problem
\begin{align*}
  {\matFont{R}}_k = \argmin_{{\matFont{R}}_k \in\text{SO(3)}} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2 + \lambda a_k\| {\matFont{R}}_k \hat{{\vecFont{n}}}_k - {\vecFont{t}}_k \|_2^2
\end{align*}
The above optimization is an instance of the \emph{orthogonal Procrustes} which finds the best rotation matrix ${\matFont{R}}_k$ to map a set of vectors (${\vecFont{e}}_{ij}, \hat{{\vecFont{n}}}_k$) to another set of vectors (${{\vecFont{e}}'}_{ij}, {\vecFont{t}}_k$). We can re-write it into a more compact expression as:
\begin{align}\label{equ:localStep}
  &{\matFont{R}}_k^\star = \argmax_{{\matFont{R}}_k \in \text{SO}(3)} \ \Tr({\matFont{R}}_k {\matFont{X}}_k) \\
  &{\matFont{X}}_k =
  \begin{bmatrix}
    {\matFont{E}}_k & \hat{{\vecFont{n}}}_k
  \end{bmatrix}
  \begin{bmatrix}
    {\matFont{W}}_k & \vspace{3.5pt}\\
    & \lambda a_k
  \end{bmatrix}
  \begin{bmatrix}
    {{\matFont{E}}'}_k^\top  \vspace{2pt}\\
    {\vecFont{t}}_k^\top
  \end{bmatrix}.
\end{align}
where ${\matFont{W}}_k$ is a $|\mathcal{N}_k|$-by-$|\mathcal{N}_k|$ diagonal matrix of the cotangent weights $w_{ij}$, ${\matFont{E}}_k$ and ${{\matFont{E}}'}_k$ are $3$-by-$|\mathcal{N}_k|$ matrices concatenating the edge vectors of the face one-ring at the rest and deformed states, respectively.
One can then derive the optimal ${\matFont{R}}_k$ from the SVD of ${\matFont{X}}_k = \mathcal{U}_k \sum_k \mathcal{V}_k^\top$
\begin{align}\label{equ:optRotation}
  {\matFont{R}}_k = \mathcal{V}_k \mathcal{U}_k^\top,
\end{align}
up to changing the sign of the column of $\mathcal{U}_k$ so that $\det({\matFont{R}}_k ) > 0$.

\section{Global Step with $E_{\textsc{arap}}\xspace$} \label{app:globalStep}
The global step updates the deformed vertex positions $\widetilde{{\matFont{V}}}$ from a fixed set of rotations ${\matFont{R}}$ obtained via the local step. This boils down to solving the following problem
\begin{align*}
  \widetilde{{\matFont{V}}}^\star = \argmin_{\widetilde{{\matFont{V}}}} \sum_{k \in V} &\sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2
\end{align*}
We can expand this energy as
\begin{align*}
  &\sum_{k \in V} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2 \\
  &\quad = \sum_{k \in V} \sum_{i,j \in f_k} w_{ij} {{\vecFont{e}}'}_{ij}^\top {{\vecFont{e}}'}_{ij} - 2 w_{ij} {{\vecFont{e}}'}_{ij}^\top {\matFont{R}}_k {\vecFont{e}}_{ij} + \text{constant}
\end{align*}
It is often convenient to express the summation in terms of matrices.
We introduce a directed incidence matrix ${\matFont{A}}_k$ with size $|V|$-by-$|\mathcal{N}_k|$ to  represent the edge vectors in $\mathcal{N}_k$ as ${\matFont{V}}^\top {\matFont{A}}_k$, and we use ${\matFont{M}}_k$ to represent a $|\mathcal{N}_k|$-by-$|\mathcal{N}_k|$ diagonal matrix of the weights $w_{ij}$.
Then we can re-write the energy in terms of matrices as
\begin{align}\label{equ:arapGlobalStep}
  &\sum_{k \in V} \Tr({\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}} \widetilde{{\matFont{V}}}^\top {\matFont{A}}_k) - 2 \Tr({\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}} {\matFont{R}}_k {\matFont{V}}^\top {\matFont{A}}_k ) \nonumber\\
  &= \sum_{k \in V} \Tr(\widetilde{{\matFont{V}}}^\top {\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}}) - 2 \Tr({\matFont{R}}_k{\matFont{V}}^\top{\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}}) \nonumber\\
  &= \Tr\Big(\widetilde{{\matFont{V}}}^\top \Big(\textstyle\sum_k {\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top \Big) \widetilde{{\matFont{V}}} \Big) -2 \Tr\Big(\Big(\textstyle\sum_k {\matFont{R}}_k{\matFont{V}}^\top{\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top\Big) \widetilde{{\matFont{V}}}\Big) \nonumber\\
  &= \Tr(\widetilde{{\matFont{V}}}^\top {\matFont{Q}} \widetilde{{\matFont{V}}}) -2 \Tr({\matFont{R}} {\matFont{K}} \widetilde{{\matFont{V}}}),
\end{align}

where ${\matFont{R}} = \{ {\matFont{R}}_k \}$ is the concatenation of all the rotations, ${\matFont{Q}}$ is a $|V|$-by-$|V|$ symmetric matrix, and ${\matFont{K}}$ is a $|9V|$-by-$|3V|$ matrix stacking the constant terms which can be computed during the precomputation. We can then find the optimal $\widetilde{{\matFont{V}}}$ by solving a linear system
\begin{align*}
  {\matFont{Q}} \widetilde{{\matFont{V}}} = {\matFont{K}}^\top {\matFont{R}}^\top
\end{align*}
As we know from \cite{SorkineA07}, ${\matFont{Q}}$ is the cotangent Laplacian \cite{PinkallP93}. We can pre-factorize ${\matFont{Q}}$ to speed up runtime performance. With these pieces in hand, we can minimize our energy \refequ{arapNormalDriven} by iteratively performing the local and the global steps (see \refalg{normalDriven}).

\section{Generalize to $E_{\textsc{farap}}\xspace$ and $E_{\textsc{acap}}\xspace$}\label{app:regularizations}
Changing the regularization from $E_{\textsc{arap}}\xspace$ to the membrane-only regularization $E_{\textsc{farap}}\xspace$ (\refequ{farap}) requires to re-define ${\matFont{R}}$ on each face and change the set of edge vectors to the three edge vectors of a triangle. These changes would lead us to replace the ${\matFont{E}}_k, {{\matFont{E}}'}_k$ in the local step \refequ{localStep} to the three edge vectors of a face, and $a_k$ to the face area. In the global step, one only needs to update the incidence matrix ${\matFont{A}}_k$ in \refequ{arapGlobalStep} to a $|V|$-by-$|f_k|$ matrix containing the three edge vectors information.

Deploying the as-conformal-as-possible regularization $E_{\textsc{acap}}\xspace$ (\refequ{acap}) changes the local step to solve an instance of the isotropic orthogonal Procrustes problem, where an analytical solution has been derived in \cite{schonemann1970fitting}. In short, one can obtain the optimal rotation the same way as \refequ{optRotation}, and compute the optimal scaling $s_k$ analytically as
\begin{align*}
  s_k = \frac{\Tr({\matFont{W}}_k {{\matFont{E}}'}^\top_k {\matFont{R}}_k {\matFont{E}}_k) + \lambda a_k \hat{{\vecFont{n}}}_k^\top {\vecFont{t}}_k}{\Tr({\matFont{W}}_k {\matFont{E}}^\top_k {\matFont{E}}_k) + \lambda a_k \hat{{\vecFont{n}}}_k^\top {\vecFont{t}}_k}.
\end{align*}
When assembling the matrices for the global step, using $E_{\textsc{acap}}\xspace$ would require replacing ${\matFont{R}}_k$ with $s_k{\matFont{R}}_k$.

\section{Projective Dynamics for Dynamic Target Normals}\label{app:PD}
We draw inspiration from \emph{projective dynamics} \cite{BouazizMLKP14} to handle cases where the target normal $T$ is a function of output geometry $B'$. Let us first define
\begin{align*}
  T = \argmin E_N(\widetilde{{\matFont{V}}})
\end{align*}
as a minimizer of an energy $E_N$ defined on the output shape. In our cases, $E_N$ could be the distance to the closet normals or the developable energy \cite{SteinGC18}. With this definition, we re-write \refequ{arapNormalDriven} as
\begin{align*}
  &\min_{\widetilde{{\matFont{V}}},{\matFont{R}}}\ \sum_{k \in V} \sum_{i,j \in\mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2  + \lambda a_k \| {\matFont{R}}_k \hat{{\vecFont{n}}}_k  - {\vecFont{t}}_k \|_2^2,\\
  &\text{ subject to}\ T = \argmin E_N(\widetilde{{\matFont{V}}})
\end{align*}
This reformulation allows us to directly deploy the projective dynamics solver by first projecting $T = \{ {\vecFont{t}}_k\}$ to the ``constraint'' $E_N$, fixing ${\vecFont{t}}_k$, and solving the original problem as \refequ{arapNormalDriven} via the local/global solver to get $\widetilde{{\matFont{V}}}$ at the next iteration. We then iterate this procedure (see \refalg{normalDriven}) until convergence.
This expression enables us to plug-and-play different $E_N$ for different modeling objectives.

\section{Normal Driven Developable Surfaces} \label{app:developable}
Our normal-driven editing can be used to create developable surfaces by specifying a set of target normals that are developable.
Stein et al. \cite{SteinGC18} propose a characterization of discrete developability based on face normals of a vertex one-ring.
In short, if all the one-ring face normals correspond to a common plane or two planes, then this local one-ring is piecewise developable.

With this characterization, we can easily get a set of ``developable'' face normals by (1) visiting all the one-ring faces of a vertex, (2) performing a small principle component analysis on the face normals for each one-ring, and (3) projecting the normals to one or two common planes by zeroing out the components correspond to the smallest eigenvalues.
By using a different threshold to decide whether to zero out the smallest or the smallest two components, we can control the amount of creases in the developable approximation (see \reffig{control_developable}).
As each face will receive three (possibly) different developable normals from the previous procedure, we simply average them to get the target face normals.
We perform this developable normal computation at each iteration in parallel, which corresponds to the Line 9 of \refalg{normalDriven}.