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 [1]. 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 [1]. The fact that second term in each RHS is a unit complex number plays a simplifying role. The results are in the next section.
Throughout this section we assume that \(L
\in [0,2\pi]\). With lots of help from Macaulay2 [2], 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:
|
|
|
|
|
|
|
|
|
|
|
|
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.
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 [3]. 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.
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().
R version 4.5.0 (2025-04-11 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows 11 x64 (build 26100) Matrix products: default LAPACK version 3.12.1 locale: [1] LC_COLLATE=C [2] LC_CTYPE=English_United States.utf8 [3] LC_MONETARY=English_United States.utf8 [4] LC_NUMERIC=C [5] LC_TIME=English_United States.utf8 time zone: America/Los_Angeles tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] flextable_0.9.7 polarzonoid_0.1-2 loaded via a namespace (and not attached): [1] zip_2.3.3 cli_3.6.5 knitr_1.50 [4] rlang_1.1.6 xfun_0.52 textshaping_1.0.1 [7] gdtools_0.4.2 jsonlite_2.0.0 data.table_1.17.2 [10] glue_1.8.0 xslt_1.5.1 openssl_2.3.2 [13] askpass_1.2.1 V8_6.0.3 htmltools_0.5.8.1 [16] ragg_1.4.0 sass_0.4.10 fontquiver_0.2.1 [19] rmarkdown_2.29 grid_4.5.0 evaluate_1.0.3 [22] jquerylib_0.1.4 fastmap_1.2.0 yaml_2.3.10 [25] lifecycle_1.0.4 compiler_4.5.0 officer_0.6.8 [28] fontLiberation_0.1.0 Rcpp_1.0.14 katex_1.5.0 [31] systemfonts_1.2.3 digest_0.6.37 R6_2.6.1 [34] curl_6.2.2 bslib_0.9.0 logger_0.4.0 [37] uuid_1.2-1 fontBitstreamVera_0.1.1 tools_4.5.0 [40] equatags_0.2.1 xml2_1.3.8 cachem_1.1.0