The focus of this vignette is the zonohedron() constructor and specifically its tolerance argument e2, whose default value is 1.e-10.
One goal of the zonohedra package is to handle all possible zonogon facets, not just the parallelograms in the generic case. The input to the constructor is matrix whose columns are the generators of the zonohedron. The generators of a specific facet span a plane, and adding another generator increases the span to all of \(\mathbb{R}^3\). Stated another way, the set of generators of a specific facet has rank 2, and is maximal with respect to this property. So a naive way of determining the facets is to examine all subsets of the generators and determine whether each one has this property.
This is hopelessly impractical. Moreover, although the rank function is well-defined for matrices with numbers in \(\mathbb{R}\), it is not computationally meaningful for floating-point numbers. For example, if a set of floating-point vectors spans the xy-plane, their rank is unambiguously 2; the smallest singular value is 0. But if the set is given a random rotation, the smallest singular value will be very small, but non-zero. Some sort of tolerance is needed.
The central dogma is that there are vector generators in \(\mathbb{R}^3\) that are very close to the given (dyadic rational floating point) vectors, and have actual rank 2. The package does a feasibility test that the floating point generators could have come from true real vectors. This test comes from the axioms of matroid theory.
The facet-finding method chosen for zonohedron() does not use rank, but it also requires a tolerance - the argument e2.
zonohedron() are:
e0 is used here
e1 is used here. Every set of two distinct generators \(\{ v_i, v_j \}\) now has rank 2, so their cross-product \(v_i \times v_j \neq 0\).
e2 as a “pseudo-angular” threshold. Special measures are taken so that vector \(u_{i,j}\) is considered identical to \(-u_{i,j}\).
If \(E\) is changed to be just a set of abstract points, then an integer-valued function defined on subsets of \(E\) that satisfies the axioms (R1), (R2), and (R3) defines a matroid on the ground set \(E\). The rank of the matroid is defined to be \(r(E)\). We mostly follow references [1] and [2].
A given matroid \(M\) may not be represented by a set of vectors in \(\mathbb{R}^n\). But if it is, we say that \(M\) is representable over \(\mathbb{R}\). We also say that \(M\) is a vector matroid.
From (R1) it follows that a point has rank 0 or 1. A point of rank 0 is called a loop; in a vector matroid a loop corresponds to the 0 vector. A multiple group is a subset of size 2 or more, which has rank 1, and with all points of rank 1, and which is maximal. In a vector matroid a multiple group is a maximal set of 2 or more non-zero vectors that are all multiples of each other.
A simple matroid is a matroid with no loops or multiple groups.
A rank function is defined for every subset of \(E\), and is much too large to deal with directly. Matroid theory provides more efficient alternatives.
In a matroid \(M\) on a ground set \(E\), a hyperplane is a maximal subset \(H \subseteq E\) with \(r(H)=r(E)-1\).
One can show that the set of hyperplanes has these properties:For a proof see [1] p. 39.
Conversely, if a collection of subsets of \(E\) satisfies the axioms (H1), (H2) and (H3), then the collection defines a valid rank function and a matroid on \(E\). To do this, first define the corank function \(c()\) by: \[\begin{equation} c(A) := \max \Bigl\{ k : \text{there are hyperplanes } H_1,..., H_k \text{ where for all } j, A \subseteq H_j \text{ and } H_1 \cap ... \cap H_{j-1} \nsubseteq H_j \Bigr\} \end{equation}\] And now define \(r(A) := c(\varnothing) - c(A)\). This function \(r()\) satisfies the axioms (R1), (R2), and (R3). The above formula appears in [2] p. 306, without a proof.
Given a collection of hyperplanes, checking the hyperplane axioms (H1), (H2), and (H3) is more efficient than checking the rank function axioms (R1), (R2), and (R3), but still too time-consuming in practice.
In a matroid \(M\) on a ground set \(E\), a circuit is a subset \(C \subseteq E\) with \(r(C)=|C|-1\) and \(r(C - x) = r(C)\) for all \(x \in C\).
One can show that the set of circuits has these properties:For a proof see [1] p. 9.
Conversely, if a collection of subsets of \(E\) satisfies the axioms (C1), (C2) and (C3), then the collection defines a valid rank function and a matroid on \(E\). \[\begin{equation} r(A) := |A| - \max \Bigl\{ k : \text{there are circuits } C_1,..., C_k \text{ where for all } j, C_j \subseteq A \text{ and } C_j \nsubseteq C_1 \cup ... \cup C_{j-1} \Bigr\} \end{equation}\] This formula appears in [2] p. 306, without a proof.
A circuit of size 1 is a loop. A circuit of size 2 is a pair of points in a multiple group. Recall that simple matroid is a matroid with no loops or multiple groups. Thus, a simple matroid is a matroid with no circuits of size 1 or 2.
In this section we derive an efficient way to check the hyperplane axioms, but only in the case when the matroid rank is 3.
Given an integer \(d \ge 1\) a \(d\)-partition of \(E\) is a collection of subsets of \(E\), called blocks, with these properties:One can show that the blocks of a \(d\)-partition satisfy the hyperplane axioms (H1), (H2), and (H3). For a proof see [1] p. 40. The resulting matroid on \(E\) is called a paving matroid and has rank \(d{+}1\). Note that the 3 properties of a \(d\)-partition can be checked efficiently.
Theorem A matroid of rank \(r \ge 2\) is a paving matroid if and only if every circuit has size \(r\) or greater.
Proof See [1], p. 40.
Theorem A simple matroid \(M\) of rank 3 is a paving matroid.
Proof (trivial) Since \(M\) is simple no circuit has size 1 or 2. Therefore every circuit has size 3 or greater. By the previous theorem, \(M\) is paving. \(\square\)
Given a set of proposed hyperplanes for a matroid of rank 3, we finally have an efficient way to check the hyperplane axioms, by checking the \(d\)-partition block axioms instead.For the hyperplane simplification in item 1, the number of hyperplanes is preserved, but all loops are removed, and every generator except one from each multiple group are removed.
To summarize, let \(E\) be a finite set of floating point 3D vectors, with no vector equal to 0 and no vector a multiple of another (with tolerances). The vectors generate a zonohedron. A collection of subsets of \(E\) is then computed, with each subset coplanar, or very close to coplanar using the tolerance parameter e2 discussed above. Each subset is the proposed set of generators of a facet of the generated zonohedron, and all facets are represented. These subsets are proposed as the hyperplanes of a matroid. We have shown that:
If \(E\) can be slightly perturbed to a set of actual real vectors \(E' \subset \mathbb{R}^3\), so that the rank of each real hyperplane is 2, and is maximal w.r.t. this property, then these hyperplanes satisfy properties (D1), (D2), and (D3).
In the software package, we use the contrapositive form:
If these proposed hyperplanes do not satisfy (D1), (D2), and (D3), then the hyperplanes do not form a valid matroid, and \(E\) cannot be slightly perturbed to satisfy the desired rank=2 property.
Even if the matroid is valid, the perturbation \(E'\) may not exist, because the matroid might not be representable over the real numbers \(\mathbb{R}\). A classical example is the Fano plane matroid on 7 points with 7 hyperplanes. It has just too many hyperplanes, see [3].
Nevertheless, we conjecture that such non-representable matroids cannot occur in practice.
Conjecture If the hyperplanes for the floating point set \(E\) are computed following the procedure in the Introduction, and the tolerance e2 (depending on \(E\)) is sufficiently small, then a perturbation \(E' \subset \mathbb{R}^3\) representing the matroid exists.
This statement is theoretical in nature, since real numbers in \(\mathbb{R}\) cannot be represented exactly.
The conjecture is true in some simple cases. Before exploring this, call the hyperplanes of size 2 the trivial hyperplanes. Note that for the Fano plane matroid, all 7 hyperplanes are size 3 and non-trivial.
Suppose that all the hyperplanes are trivial, so the matroid is uniform and all the facets of the zonohedron are parallelograms. Then no perturbation is needed at all; the given vectors (with dyadic rational numbers) already represent. This is the case for 7 of the 13 classical zonohedra in classics.genlist. And it is also the case for the generators in colorimetry.genlist[[3]].
Now suppose that the matroid has only 1 non-trivial hyperplane. Then there are 3 or more generators that (approximately) span a plane, and all the other generators are far from the plane. Perturb this plane to the “best fit” linear plane \(P\) to these generators where \(P \subset \mathbb{R}^3\), and then project them onto \(P\). If this perturbation accidentally creates non-trivial hyperplanes with the other generators, then just perturb the other generators to get the original matroid. An example is the matroid generated by colorimetry.genlist[[2]], which has 1 non-trivial hyperplane with 50 generators.
Now suppose that all the non-trivial hyperplanes are disjoint. Then we can repeat the procedure in the previous paragraph for each hyperplane. Since the hyperplanes are disjoint, there is no “interaction” between them. An example is the matroid generated by colorimetry.genlist[[1]], which has 2 disjoint non-trivial hyperplanes with sizes 3 and 26.
Now suppose that the non-trivial hyperplanes intersect in a single generator. We can perform a “constrained best fit” perturbation for each plane \(P\), where the constraint is that that single generator is in the plane. An example is the matroid generated by classics.genlist[[5]], which has 2 non-trivial hyperplanes: \(\{1, 3, 4\}\) and \(\{2, 3, 5\}\). The generated zonohedron is the rhombo-hexagonal dodecahedron.
More simple cases can be listed by mixing the above, but we cannot find a general proof of the conjecture.
R version 4.3.0 (2023-04-21 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045) Matrix products: default 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 loaded via a namespace (and not attached): [1] digest_0.6.31 R6_2.5.1 fastmap_1.1.1 xfun_0.39 [5] cachem_1.0.8 knitr_1.42 htmltools_0.5.5 rmarkdown_2.21 [9] cli_3.6.1 sass_0.4.6 jquerylib_0.1.4 compiler_4.3.0 [13] tools_4.3.0 evaluate_0.21 bslib_0.4.2 yaml_2.3.7 [17] rlang_1.1.1 jsonlite_1.8.4