Implicitization"
In polarzonoid: Compute Maps and Properties of Polar Zonoids

```{css, echo=FALSE} body { max-width: 725px; / make wider, default is 700px / }

```r
knitr::opts_chunk$set(echo = TRUE)
old_opt = options( width=144 )

Introduction

library(polarzonoid)

In the User Guide vignette it is shown that there is are homeomorphisms \begin{equation} A_n ~~ \longleftrightarrow ~~ \partial Z_n ~~ \longleftrightarrow ~~ \mathbb{S}^{2n} \end{equation} where $A_n$ is the space of $n$ or fewer pairwise disjoint arcs in the circle, and $Z_n$ is the polar zonoid in $\mathbb{R}^{2n+1}$. The forward direction $\partial Z_n ~~ \to ~~ \mathbb{S}^{2n}$ is the trivial projection onto the unit sphere centered at $(0,...0,\pi)$. But the reverse direction requires intersecting a ray with $\partial Z_n$, which is non-trivial. The function in the package that performs this intersection is boundaryfromsphere(). It expresses $\partial Z_n$ as the level set of a real-valued function, and uses uniroot() to solve for the intersection point.

In this vignette, we outline how $\partial Z_n$ is defined parametrically and how the real-valued function is derived. This process is called implicitization, see @cox2008. In the current version of the package, the implicitization has only been computed for $n$=1,2, and 3.

Throughout this vignette we think of $\partial Z_n \subseteq \mathbb{R}^{2n+1}$ and identify the latter with points $(z_1,...,z_n,L)$ where $z_i$ are complex numbers and $L$ is real.

Assume first that $n{=}1$ so there is a single arc in the circle, i.e. a point in $A_1$. We parameterize the arc by its first and last endpoints $u_1$ and $v_1$, both on the unit circle in the complex plane, and proceeding counterclockwise from $u_1$ to $v_1$. Let the length of the arc be $l_1$, so that $v_1 = u_1 e^{i l_1}$. It is not hard to show that the map $A_1 \to \partial Z_1$ is given by \begin{equation} (u_1,v_1) \mapsto ( ~(v_1{-}u_1)/i, ~ l_1 ~) ~=~ (z_1,L) \end{equation} where $z_1 := (v_1{-}u_1)/i$. Similarly, when $n{=}2$, it is easy to show that the map $A_2 \to \partial Z_2$ is given by: \begin{equation} (u_1,v_1,u_2,v_2) \mapsto ( ~ \left( v_1{-}u_1 + v_2{-}u_2 \right)/i, ~ \left( v_1^2{-}u_1^2 + v_2^2{-}u_2^2 \right)/(2i), ~ l_1{+}l_2 ~) ~=~ (z_1,z_2,L) \end{equation} and so on for $n \ge 3$. Since $u_i$ and $v_i$ are unit vectors, $u_i \overline{u_i} = 1$ and $v_i \overline{v_i} = 1$ where the overline denotes complex conjugate. It follows easily that: \begin{align} \overline{z_1} &= \big( u_2 v_2 (v_1 - u_1) + u_1 v_1(v_2 - u_2) \big) \left( -i / (u_1 u_2 v_1 v_2) \right) \ \overline{z_2} &= \big( u_2^2 v_2^2 (v_1^2 - u_1^2) + u_1^2 v_1^2(v_2^2 - u_2^2) \big) \left( -i / (2 u_1^2 u_2^2 v_1^2 v_2^2) \right) \end{align}

And using Euler's formula, one can show: \begin{align} \sin(L/2) &= \frac{v_1 v_2 - u_1 u_2}{2i} ~~ \frac{e ^{-iL/2}}{u_1 u_2} \ \cos(L/2) &= \frac{v_1 v_2 + u_1 u_2}{2} ~~ \frac{e ^{-iL/2}}{u_1 u_2} \end{align}

The expressions for $n{=}1$ and $n{=}3$ are similar. These right-hand-sides are "rational enough" that algebraic geometry can be used for implicitization. We use the Rational Implicitization algorithm in Chapter 3, Section 3 of @cox2008. The fact that second term in each RHS is a unit complex number plays a simplifying role. The results are in the next section.

The Equations for $n$ = 1,2, and 3 arcs

Throughout this section we assume that $L \in [0,2\pi]$. With lots of help from Macaulay2 @M2, followed by manual rearrangement, one can show that for $n$ = 1,2, and 3 that \begin{equation} (z_1,...,z_n,L) \in \partial Z_n ~~~ \text{ iff } ~~~ | \gamma_n(z_1,...,z_n,L) | = \rho_n(z_1,...,z_n,L)
\end{equation} and \begin{equation} (z_1,...,z_n,L) \in Z_n ~~~ \text{ iff } ~~~ | \gamma_n(z_1,...,z_n,L) | \le \rho_n(z_1,...,z_n,L)
\end{equation} where $\gamma_n()$ is a complex-valued function, and $\rho_n()$ is a real-valued function. The details of this are complicated, especially for $n{=}3$, and are omitted. For $n{=}4$, Macaulay2 runs out of memory and cannot finish; but I conjecture that these two functions exist for all $n$. For small $n$, they are given in this table:

library( flextable )
requireNamespace('equatags',quietly=TRUE)

rhs <- c(
  "\\gamma_n(z_1,...,z_n,L)",
  "z_1",
  "\\cos(L/2)z_1^2 ~-~ 2 \\sin(L/2)z_2",
  "12 ( 4\\sin^2(L/2) ~-~ |z_1|^2 ) z_3 \\newline { ~+~ (8\\cos^2(L/2) - |z_1|^2 + 4) z_1^3 } \\newline ~-~ 48 \\sin(L/2)\\cos(L/2) z_1 z_2 ~+~ 12 \\overline{z_1} z_2^2"
  )

lhs = c(
    "\\rho_n(z_1,...,z_n,L)",
    "2\\sin(L/2)", 
    "4\\sin^2(L/2) ~-~ |z_1|^2", 
    "{ 6\\sin(L/2) (|z_1|^4 - 8|z_1|^2 - 4|z_2|^2 ) }  \\newline { ~+~ 12\\cos(L/2)(\\overline{z_1}^2z_2 + z_1^2\\overline{z_2}) }   \\newline ~+~ 96\\sin^3(L/2)" )

df <- data.frame( n=c('n',1:3), gamma=rhs, rho=lhs ) # ; df

ft <- flextable(df)
ft <- compose( x=ft, j="n", value = as_paragraph(as_equation(n, width=0.5)))
ft <- compose( x=ft, j="gamma", value = as_paragraph(as_equation(gamma, width=4.2)))
ft <- compose( x=ft, j="rho", value = as_paragraph(as_equation(rho, width=4.2)))
ft <- align( ft, align = "center", part = "all")
ft <- width( ft, j=c(2,3), width=4.2 )
ft <- delete_part( ft, part="header" )

#  add borders
thick_border = fp_border_default(color="black", width = 2)
thin_border = fp_border_default(color="gray", width = 1)
ft <- border_inner_h( ft, part="all", border = thin_border )
ft <- border_inner_v( ft, part="all", border = thin_border )
ft <- border_outer( ft, part="all", border = thick_border )
ft <- hline( ft, i=1, border=thick_border )
ft <- vline( ft, i=1, j=c(1,2), border=thick_border )

ft

The functions $\rho_1()$ and $\rho_2()$ are obviously real, and $\rho_3()$ is real because $\overline{z_1}^2z_2 + z_1^2\overline{z_2}$ is real.

The case $n{=}1$ is in $\mathbb{R}^3$ and easy to visualize: \begin{equation} (z_1,L) \in \partial Z_1 ~~~ \text{ iff } ~~~ |z_1| = 2\sin(L/2) \end{equation} This is a surface of revolution - a sine wave around the L axis. It is like a football with poles at (0,0) and (0,$2\pi$). Note that these 2 poles are exactly the points where the right side of the equality vanishes.

If we define $f(z_1,...,z_n,L) := | \gamma_n(z_1,...,z_n,L) | - \rho_n(z_1,...,z_n,L)$ then $f()$ vanishes on $\partial Z_n$, is negative in the interior of $Z_n$, and is positive outside $Z_n$. This level-set function $f()$ is used in boundaryfromsphere() to compute the intersection of a ray based at $(0,\pi)$ and $\partial Z_n$. The gradient of $f()$ is used in arcsfromboundary(); see the next section.

Supporting Hyperplanes of $\partial Z_n$

For $x = (z_1,...,z_n,L) \in \partial Z_n$ evaluating the function $\partial Z_n \to A_n$ requires the outward normal of a supporting hyperplane to $\partial Z_n$ at $x$. The coordinates of this normal become the coefficients of a trigonometric polynomial, whose roots are the endpoints of arcs. If $\partial Z_n$ is smooth at $x$ the supporting hyperplane is unique, and is in fact the gradient of the level-set function $f()$ at $x$. Also in this case the trigonometric polynomial has the full set of $2n$ roots, defining $n$ arcs.

But unfortunately, $\partial Z_n$ is NOT smooth at all $x$. We saw this already with $n{=}1$, where the surface is not smooth at the 2 poles. At the poles, there are infinitely many supporting planes, the trigonometric polynomial has 0 roots, the "south pole" maps to the empty arc, and the "north pole" to the full circle.

The function $\partial Z_n \to A_n$ is implemented in arcsfromboundary(). This R function has to deal with the non-smooth case. It turns out that $\partial Z_n$ is smooth at $x$ iff $\rho_n(x) > 0$. We saw this already with $n{=}1$. So starting at $k{=}1$ and incrementing up to $n$, it tests $\rho_k(z_1,...,z_k,L)$. If $\rho_k(z_1,...,z_k,L) < \epsilon$, it returns a collection of $k{-}1$ arcs in $A_n$. Otherwise the surface is smooth at $x$ and it returns $n$ arcs. Recall that $A_n$ is the space of $n$ or fewer arcs. arcsfromboundary() first calculates the outward pointing normal at p. The the gradient of the level-set function $f()$ is used to compute this normal; the expressions for the gradient are more complicated than the ones in the above table, and are omitted here. This coordinates of this normal are the coefficients of a trigonometric polynomial whose roots are calculated with trigpolyroot(). These roots are the endpoints of the arcs.

NOTE: The differences $X_i := \partial Z_i - \partial Z_{i-1}$ almost surely form a Thom-Mather stratification of $\partial Z_n \subset \mathbb{R}^{2n+1}$, but I have not checked this. The conditions for such a stratification are complex, and are in @ThomMather. Each $X_i$ is the image of exactly $i$ arcs. I think that the above functions $\rho_i()$, after restriction to $X_i$, can serve as the distance functions in the stratification.

Summary

The forward functions: \begin{equation} A_n ~~ \longrightarrow ~~ \partial Z_n ~~ \longrightarrow ~~ \mathbb{S}^{2n} \end{equation} are straightforward. In the package, they are implemented as boundaryfromarcs() and spherefromboundary() respectively.

But the inverse functions: \begin{equation} A_n ~~ \longleftarrow ~~ \partial Z_n ~~ \longleftarrow ~~ \mathbb{S}^{2n} \end{equation} use convergent numerical iteration. The first one, boundaryfromsphere(), uses the implicit function defined above, and this version of the package requires that $n$ = 1, 2, or 3. The second one, arcsfromsphere(), solves for the roots of a trigonometric polynomial and uses base::polyroot().