Epsilon local rigidity and numerical algebraic geometry
EEpsilon local rigidity and numerical algebraic geometry
Andrew Frohmader and Alexander Heaton University of Wisconsin, Milwaukee Max Planck Institute for Mathematics in the Sciences, Leipzig, and TechnischeUniversit¨at BerlinFebruary 12, 2020
Abstract
A well-known combinatorial algorithm can decide generic rigidity in the plane by determining if thegraph is of Pollaczek-Geiringer-Laman type. Methods from matroid theory have been used to prove otherinteresting results, again under the assumption of generic configurations. However, configurations arisingin applications may not be generic. We present Theorem 5 and its corresponding Algorithm 1 whichdecide if a configuration is ε -locally rigid, a notion we define. A configuration which is ε -locally rigidmay be locally rigid or flexible, but any continuous deformations remain within a sphere of radius ε inconfiguration space. Deciding ε -local rigidity is possible for configurations which are smooth or singular,generic or non-generic. We also present Algorithms 2 and 3 which use numerical algebraic geometry tocompute a discrete-time sample of a continuous flex, providing useful visual information for the scientist. Keywords: Numerical algebraic geometry, Real algebraic geometry, Rigidity, Kinematics, Mechanism Mobil-ity, Homotopy Continuation. 2010 Mathematics Subject Classification: 70B15, 65D17, 14Q99.
Consider a graph with n nodes labelled by [ n ] = { , , . . . , n } , and with edge set E ⊂ (cid:0) [ n ]2 (cid:1) . A map calledthe initial configuration p : [ n ] → R d embeds this graph in some R d . We make precise definitions in Section2, but the basic idea is to consider the edge distances (cid:96) ij between nodes connected by an edge, but also thedistances (cid:99) (cid:96) ij between nodes that are not connected by an edge. If the nodes can move so that the distances (cid:96) ij for { i, j } ∈ E remain constant, but some distance (cid:99) (cid:96) ij for { i, j } / ∈ E changes, then we say the configuration p is flexible . If no such continuous motion exists, we say the graph is locally rigid . We can always translateor rotate the graph inside R d , called a rigid motion, which keeps all pairwise distances constant. Thus, todecide local rigidity, we need to establish if there are deformations of the embedded graph besides the rigidmotions .It is often difficult to determine conclusively whether a given configuration p admits deformations whichpreserve edge-lengths, yet are not rigid motions (see [1] where they show this problem is coNP-hard). Suchdeformations preserving the (cid:96) ij but changing some non-edge length (cid:99) (cid:96) ij are called flexes . We demonstratethese ideas with an historically important example (appearing in the first 3 competing patents for tensegrity[21]) which we call the 3-prism, whose rigidity is well-known. Any generic configuration of the 3-prism isinfinitesimally rigid, and therefore locally rigid (see the well-known Theorem 1 below). But there are manynon-generic configurations which lie on the singular locus of the polynomial system of member constraints.The left side of Figure 1 illustrates a sequence of configurations nearby such a singular p appearing to1 a r X i v : . [ m a t h . M G ] O c t eform the 3-prism. These configurations were computed numerically using Algorithm 3 below. Indeed, upto numerical tolerance, all the distances (cid:96) ij between nodes connected by an edge are the same, yet severalnon-edge distances (cid:99) (cid:96) ij change. Theorem 1 fails to apply since p is singular and infinitesimal flexes existFigure 1: Numerical flex of the 3-prism(right side of Figure 1). This infinitesimal flex was used as input to our Algorithm 3. Given only thisnumerical evidence we may be tempted to think that p is flexible . However, using our Theorem 5 andAlgorithm 1, we can conclusively determine that the numerical configurations illustrated in Figure 1 are notpart of a legitimate flex. To be clear, our results do not prove local rigidity, but they prove something thatmay be good enough in practice, ε -local rigidity. In contrast to many methods and rules applied in rigiditytheory, the methods presented here apply equally well to any configuration, be it exceptional or generic,singular or smooth.Hence, the main results of this article are Theorem 5 and its associated Algorithm 1. These concern ε -localrigidity, which we depict in Figure 2. We make precise definitions in Section 4, but briefly, a configuration p is ε -locally rigid if we can certify that any continuous flexes are extremely small, and cannot go far . Namely,we will certify that any continuous flex through p stays within an ε -ball about p within configuration space R nd . Since the choice of ε > ε -local rigidity is practically as good as knowinglocal rigidity. Our methods use polynomial homotopy continuation and numerical algebraic geometry [5, 27],taking advantage of theorems in real algebraic geometry [3, 13, 24, 25] which can guarantee finding real solu-tions to systems of polynomial equations, should they exist. In this way, our methods deal directly with thereal algebraic set, applying equally well to both smooth and singular configurations. In many cases ε -localrigidity may be more relevant than local rigidity, since allowing some small movement may be acceptablefor certain applications. Finally, our methods also imply Algorithms 2 and 3 which produce animations of aflex, should it exist, yielding easily-understandable visual information for the scientist. Related literature:
We briefly discuss some research related to this article. Since the simplest modelsare often the most useful, variations on this theme have been well-studied and therefore have many names: configuration, truss, bar-and-joint framework, tensegrity, linkage, assembly mode, structure, mechanism,mobility, degrees of freedom , and more. The combinatorics community has studied this problem by assumingcertain generic conditions on p . This allows for theorems that use the graph structure alone (for overviewssee [8, 31] or Chapter 61 of [12]), but requires generic assumptions. In contrast, the results of this paperapply to any configuration. The engineering and kinematics communities have also studied variations of thisproblem. In [11], a critical review of methods of mobility analysis is presented. The author enumerates 352igure 2:different approaches to calculating the mobility of a given configuration p , describing the limitations andoutright failure of the methods in various cases. Again, in contrast, the results of this article apply to anyconfiguration without assumptions. However, we cannot prove local rigidity, only ε -local rigidity.The 2019 paper [19] defines the almost rigidity of frameworks . This paper extracts information from thesingular value decomposition of the rigidity matrix for a given configuration p , using existing tests for pre-stress stability [9, 10] and semidefinite programming. Given a specific configuration p , if a certain list ofconditions is met, they provide a specific radius η > p will remain. The size of η depends on p . The key difference between their result and ours is thatour methods allow a freely chosen ε which can be adjusted according to the application. In contrast, themethods of [19] output a single radius η computed from the SVD of the rigidity matrix for p . This section collects the basic definitions and results from rigidity theory and numerical algebraic geometryneeded for this article. For general references on combinatorial rigidity theory see [8, 31] or Chapter 61 of[12]. For numerical algebraic geometry see [5, 27]. We consider a connected graph (
V, E ) with n nodes and m edges where V = [ n ] = { , , . . . , n } and an edge is written either { i, j } ∈ E , or briefly ij ∈ E , where i ∈ [ n ] and j ∈ [ n ] correspond to nodes. The graph is embedded in R d by the map p : [ n ] → R d , calledthe initial configuration . We assume that the affine span of the n nodes is d dimensional (they are not allcontained in some lower-dimensional subspace). By slight abuse of notation we specify p by a list of the nd coordinates of its nodes, denoted by a tuple ( p ik ) ∈ R nd , where i ∈ [ n ] and k ∈ [ d ], so that p ik is the k thcoordinate of the i th node. A deformation of p is a continuous path p ( t ) : [0 , → R nd which recovers p at t = 0. A rigid motion is a deformation which preserves the pairwise squared distances (cid:40) d (cid:88) k =1 ( p ik ( t ) − p jk ( t )) (cid:41) ij ∈ ( [ n ]2 ) , (1)at every time t , where ij ∈ (cid:0) [ n ]2 (cid:1) runs over all pairs of nodes, not just the nodes that are connected by anedge ij ∈ E . We associate a system of polynomial equations called the member constraints to our embeddedgraph which enforce the requirement that edges must have constant length.3 efinition 1. Let g : C nd → C m be given by g ( x ) = g ( x ) g ( x )... g m ( x ) , where g l : C nd → C gives the difference in squared length of the l th edge from that of the initial configuration p . If the l th edge is ij , then g l ( x ) = d (cid:88) k =1 ( x ik − x jk ) − d (cid:88) k =1 ( p ik − p jk ) . (2)The system of member constraints associated to p is the system of polynomial equations given by g ( x ) = 0.We will also need the corresponding algebraic set V ( g ) := { x ∈ C nd : g ( x ) = 0 } and real algebraic set V R ( g ) := V ( g ) ∩ R nd .The real algebraic set V R ( g ) is the configuration space . We can now state the key definition. Definition 2. A flex of p is a deformation p ( t ) : [0 , → R nd such that g ( p ( t )) = 0 for all t ∈ [0 ,
1] andwhich is not a rigid motion. The configuration p is called locally rigid if no flex exists.The first and easiest way to decide whether a given configuration p is locally rigid is to examine thelinearization of the polynomial map g of Definition 1. The Jacobian of g , denoted dg , is an m × nd matrixof polynomials. The matrix dg is usually called the rigidity matrix. When we evaluate dg at the point p ∈ V R ( g ), the vectors in its right nullspace N ull ( dg | p ) are sometimes called infinitesimal mechanisms , forexample in [28]. The Lie algebra of the Euclidean group of rigid motions produces (cid:0) d +12 (cid:1) linearly independentinfinitesimal mechanisms, whose span we denote RM . Any remaining vectors outside the span of theinfinitesimal rigid motions are called infinitesimal flexes , giving rise to the following decomposition: N ull ( dg | p ) = RM ⊕ F, (3)where any v / ∈ RM is an infinitesimal flex. Any configuration p for which F = 0 is called infinitesimallyrigid . The good news comes in the form of the well-known theorem Theorem 1 ([2]) . Infinitesimal rigidity implies local rigidity.Thus if rank (cid:0) dg | p (cid:1) = nd − (cid:0) d +12 (cid:1) then our configuration is infinitesimally rigid , and therefore also locallyrigid . Crucially, the converse is false - a configuration may have infinitesimal flexes but still be locally rigid.Other techniques are required to determine if infinitesimal flexes are actually realizable. At most points p ∈ V ( g ) the rank of the Jacobian will be equal to the generic rank of the polynomial map g [27, Section13.4]. But at singular points p the Jacobian dg | p will drop rank, admitting more infinitesimal mechanisms.For the 3-prism, an infinitesimal flex vector v ∈ R nd with dg | p ( v ) = 0 is pictured on the right of Figure 1by arranging its nd = 18 components as 3-vectors attached to each of the n = 6 nodes.From the discussion above, we see that infinitesimal rigidity is a sufficient condition for local rigidity but isnot necessary. The complete picture comes from the local geometry of the configuration space V R ( g ) around p . A flex is nothing more than a path through configuration space starting at p which is not a rigid motion.The configuration p is locally rigid exactly when the local real dimension of V R ( g ) is (cid:0) d +12 (cid:1) , the dimensionof the Euclidean group of rigid motions. See [30] for more details. Thus, a technique for determining thelocal real dimension of an algebraic set would completely solve the problem of determining local rigidity. Nosuch technique currently exists.A partial answer to the problem of determining local real dimension comes from the local dimension testdescribed in [30]. With minor caveats, the local dimension test determines the local complex dimension4uccessfully. The real and complex local dimensions agree whenever p is a smooth point on the algebraicset. Thus, for smooth points [30] solves the local rigidity problem. However, many of the examples we areinterested in are singular configurations. For singular p like the 3-prism of Figure 1, the real local dimensionmay differ from the complex local dimension computed in [30]. The authors acknowledge this point and leaveit open for future work.To understand the local geometry of V R ( g ) around p we rely on tools from numerical algebraic geometry - specifically polynomial homotopy continuation. Given a square system of polynomial equations (squaremeans same number of variables as equations), polynomial homotopy continuation can reliably compute all isolated solutions , in contrast to other numerical methods like Newton iteration which can only reliablycompute some solutions . Numerical algebraic geometry draws on results from algebraic geometry to providesuch guarantees. For general references see any of [5, 14, 27]. Here we describe only the basics necessary tounderstand the results of this article.Given a system g : C N → C N whose solutions we would like to discover, the key is to produce a startingsystem f : C N → C N whose solutions are known, and a homotopy H ( x, t ) : C N × C → C N such that H ( x,
1) = f ( x ) and H ( x,
0) = g ( x ). The key computational step is to solve a Davidenko initial valueproblem for the system of ODE’s resulting from differentiating H ( x ( t ) , t ) = 0 with respect to time. Here, x ( t ) is the smooth path of some known solution f ( x (1)) = 0 to some unknown solution g ( x (0)) = 0. Solvingthe ODE from the initial conditions at H ( x (1) ,
1) = 0 to the endpoint [27, Chapter 10] at t = 0 using anyof the standard ODE solvers is called path-tracking [27, Section 2.3]. For excellent examples see the websiteof the julia implementation HomotopyContinuation.jl , which we use in all our examples [7].A total-degree homotopy always provides a valid starting system [27, Section 8.4.1]. In many cases thereare smaller starting systems which require tracking less paths, and hence less computation [27, Section 8.4,8.5]. We demonstrate this for the 3-prism in Example 2 later. In addition to finding isolated solutions, theconcept of witness sets can be used to understand positive-dimensional components of the solution set topolynomial systems [27, Chapter 13]. Say X is a k -dimensional irreducible component of the solution set to g : C N → C m . A witness set for X is a tuple ( g, L, W ) where g is the polynomial system, L is a linear spaceof dimension N − k , and W is a list of the finitely many isolated solutions obtained by intersecting X withthe linear space L . Once a witness set for X has been computed, one can easily sample more points from X by perturbing L to nearby linear spaces L (cid:48) , following the solutions in W . The size | W | is the degree of X .Numerical algebraic geometry is a powerful tool that we use in this article. Our goal in this section is to remove rigid motions without changing the space of flexes. Often, when buildingstructures, you attach or fasten them to the ground, or to the wall, or to another structure. In that case,the correct model will treat those nodes as fixed, and not allow them to move at all. This deletes nodevariables from the equations g and columns from the Jacobian matrix dg . All our results still apply, thoughthe change of coordinates we describe presently would not be necessary.However, if you are considering a structure which you do not intend to fasten to the ground, like NASA’sSuper Ball Bot of Figure 3, then the correct model does not have variables deleted from the system g . Wewould like handle both situations equally well. With no nodes fixed, the configuration space includes thefull (cid:0) d +12 (cid:1) dimensional space of rigid motions. These rigid motions are independent of the structure beingstudied and are not interesting. We would like to remove them so we can focus on the flexes. What followsis a convenient procedure for removing the rigid motions. This can be viewed as a quotient constructionwhere we work with flex representatives after modding out by rigid motions.Observe that given a valid configuration ( x ik ) ∈ R nd which satisfies the member constraints g , any sequenceof rigid motions applied to ( x ik ) will yield another valid configuration ( (cid:99) x ik ). We use this to choose aconvenient reference frame for our initial configuration and stay with this reference frame as we follow a flex.5igure 3: Structure with no nodes fastened to the groundAs an example, we illustrate our choice of reference frame for a configuration embedded in R . Letnodes 1,2,3 be noncollinear. Translate the reference frame until node 1 is at the origin. Next, rotate thereference frame until node 2 lies on the x -axis. Finally, rotate the reference frame about the x -axis untilnode 3 has z -coordinate zero. The resulting configuration has node 1 = (0 , , x , , x , x , x = x = x = x = x = x = 0 at all points along a path through V R ( g ). We claim these constraints remove the rigid motionsbut do not change the space of flexes. To see this, note any valid configuration ( x ik ) can be sent, via rigidmotions, to another configuration ( (cid:99) x ik ) which satisfies the above constraints. We refer to this change ofcoordinates as the moving frame. The general case is presented in the following two theorems. Theorem 2.
Let the affine span of nodes 1 , . . . , d be d − p ik to (cid:99) p ik such that for i ∈ [ d ] we have (cid:99) p ik = 0 if k ≥ i . We denote the newconfiguration (cid:98) p . Proof.
We proceed by induction on the number of nodes j . Let the dimension of ambient space be fixed at d . If j = 1, apply a translation that moves node 1 to the origin. Then its new coordinates satisfy (cid:99) p k = 0 if k ≥ j = d −
1. We find a rigid motion that fixesnodes 1 , . . . , d − d such that p dd = 0. Let H be the subspace of R d spanned by the first d − K be the subgroup of SO ( R d ) that fixes H . By induction, nodes 1 , . . . , d − H . K is isomorphic to SO ( R ) so we can select a rotation r ∈ K that fixes indices 1 , . . . , d − d such that p dd = 0 as desired. Remark 1.
We view (cid:98) p as an element of R N where N = nd − (cid:0) d +12 (cid:1) by dropping the newly zero coordinates. Definition 3.
We associate to (cid:98) p a system of equations called the moving frame member constraints denoted (cid:98) g : C N → C m . As with g , this system enforces the requirement that edges must have constant length and its definition isanalogous to that of g . We simply drop the unnecessary variables. Theorem 3.
There exists a flex p : [0 , → R nd of the initial configuration p if and only if there exists aflex (cid:100) p ( t ) : [0 , → R N of the initial configuration (cid:98) p . Proof.
Order the nodes such that nodes 1 , . . . d span d − R nd . Let π : R nd → R N be thenatural projection and ι : R N → R nd be the natural injection. Let r be the rigid motion of Theorem 2 suchthat π ( r · p ) = (cid:98) p . First, say we have a flex (cid:100) p ( t ) : [0 , → R N . Then r − · ι ( (cid:100) p ( t )) is a flex of p in R nd . Nowsuppose we have a flex p : [0 , → R nd . By Theorem 2 for every point p ( t ), there exists a rigid motion r ( t )which sends p ( t ) to another point with (cid:100) p ( t ) ik = 0 if k ≥ i and i ∈ [ d ]. The continuity of p ( t ) implies r ( t ) iscontinuous. Thus π ( r ( t ) · p ( t )) is a flex of (cid:98) p .For interested readers, we make note of the connection to the method of moving frames [23]. Understoodin this sense, we are making our calculations in a coordinate cross-section .6 Epsilon local rigidity
In this section we determine whether a given initial configuration p is ε -locally rigid. We fix N = nd − (cid:0) d +12 (cid:1) .It is tempting to define a notion of ε -locally rigid in the original configuration space R nd but since anyconfiguration can be translated by an arbitrary amount, the better notion is applied to the point afterchanging to the moving frame (cid:98) p ∈ R N as in Theorem 2. This motivates the following Definition 4.
Let p be an initial configuration and (cid:98) p be the configuration in the moving frame of Theorem2. We say that p is ε -locally rigid if every flex (cid:100) p ( t ) of (cid:98) p satisfies (cid:100) p ( t ) ∈ B ε ( (cid:98) p ) for all t ∈ [0 , B ε ( (cid:98) p )is the open ε -ball centered at (cid:98) p . Remark 2.
In other words, if p is ε -locally rigid, then any positive-dimensional connected component ofthe real algebraic set V R ( (cid:98) g ) containing (cid:98) p stays within some ball about (cid:98) p . Any flex that may exist can besafely ignored if ε is sufficiently small. For all practical purposes, it is as if p is locally rigid.We now move towards Theorem 5 and Algorithm 1 to decide whether a configuration p is ε -locally rigid. Definition 5.
Let (cid:98) g = [ (cid:98) g , . . . , (cid:99) g m ] T be the moving frame member constraints associated to (cid:98) p as in Definition3. Define the polynomial system (cid:98) g ε : C N → C (cid:98) g ε = (cid:98) g + · · · + (cid:99) g m + s ε , where s ε is defined by s ε = ε − d (cid:88) k =1 n (cid:88) i =1 ( x ik − p ik ) . The system of ε - member constraints associated to (cid:98) p is given by (cid:98) g ε ( x ) = 0. We denote the correspondingalgebraic set (cid:98) V ε := { x ∈ C N : (cid:98) g ε ( x ) = 0 } . Lemma 1.
The irreducible components of (cid:98) V ε are of dimension exactly N − Proof.
By Theorem 13.4.2 of [27], the possible dimensions of irreducible components X of an algebraic set V ( f ) for f : C N → C n are bounded between N − rank f ≤ dim X ≤ N − , where the rank of f is the dimension of the closure of its image as a map, or equivalently, the generic rank ofits Jacobian. For a single, nonzero polynomial like (cid:98) g ε : C N → C we have that N − ≤ dim X ≤ N − Assumption 1.
We collect here the following list of assumptions which refer to the homotopy H ( x, λ, t )defined in Theorem 4 below.1. Let N > k > f : R N → R N − k be a polynomial system with real coefficients, with V ⊂ V ( f ) apure k -dimensional algebraic set with witness set { f, L, W } .2. Assume that the starting solutions to H ( x, λ,
1) = 0 are finite and nonsingular.3. Assume also that the number of starting solutions is equal to the maximum number of isolated solutionsto H ( x, λ,
1) = 0 as z, γ, y, α vary over C N − k × C × C N × C N − k +1 . This will be true for a nonemptyZariski open set of C N − k × C × C N × C N − k +1 .7. Assume all the solution paths defined by H starting at t = 1 are trackable. This means that for eachstarting solution ( x ∗ , λ ∗ ) there exists a smooth map ξ : (0 , → C N × C N − k +1 with ξ (1) = ( x ∗ , λ ∗ )and for all t ∈ (0 ,
1] we have ξ ( t ) is a nonsingular solution of H ( x, λ, t ).5. Assume that each solution path converges, collecting the endpoints of all solution paths in the sets E and E = π ( E ) where π ( x, λ ) = x projects onto the x coordinates, forgetting the λ coordinates. Theorem 4 (Theorem 5 of [13]) . Suppose that the conditions in Assumption 1 hold. Let z ∈ R N − k , γ ∈ C , y ∈ R N − V R ( f ), α ∈ C N − k +1 , and H : C N × C N − k +1 × C → C N − k +1 be the homotopy defined by H ( x, λ, t ) = f ( x ) − tγzλ ( x − y ) + λ ∇ f ( x ) T + · · · + λ N − k ∇ f N − k ( x ) T α T λ − (4)where f ( x ) = [ f ( x ) , . . . , f N − k ( x )] T . Then E ∩ V ∩ R N contains a point on each connected component of V R ( f ) contained in V . Theorem 5.
Let p be an initial configuration and (cid:98) g ε be according to Definition 5 above. Taking f = (cid:98) g ε inTheorem 4, we find that if conditions two through five in Assumption 1 are met, then E ∩ (cid:98) V ε ∩ R N = ∅ implies that p is ε -locally rigid. Proof.
Take f = (cid:98) g ε and V = (cid:98) V ε in the notation of Theorem 4 above. We have N − k = 1 and by Lemma1 all irreducible components are of dimension k = N −
1. Therefore, the first condition of Assumption 1 ismet. Say E ∩ (cid:98) V ε ∩ R N = ∅ but p is not ε -locally rigid. Let (cid:98) p ( t ) be a flex such that (cid:98) p (1) / ∈ B ε ( (cid:98) p ) and let P = (cid:98) p ([0 , (cid:98) p . Then (cid:16) P ∩ B ε ( (cid:98) p ) (cid:17) ∪ (cid:16) P ∩ B ε ( (cid:98) p ) c (cid:17) is a separation of P contradicting the continuity of (cid:98) p .Theorem 5 suggests the following Algorithm 1. Algorithm 1:
Epsilon local rigidity
Input:
Initial configuration p ∈ R nd , edge set E , and choice of ε > Result:
Boolean v which is true if items 2, 4, and 5 of Assumption 1 are satisfied. Boolean u whichis true if the set E ∩ (cid:98) V ∩ R N of Theorem 5 is empty, and false otherwise. Set R which maybe empty or else contains at least one point on each connected component of V R ( (cid:98) g ε ). Apply the rigid motions of Theorem 2 to p obtaining (cid:98) p ∈ R N for N = nd − (cid:0) d +12 (cid:1) . Form the systems of equations (cid:98) g ε according to Definition 5. Calculate a witness set W for the pure N − V ( (cid:98) g ε ) ⊂ C N . Produce z, γ, y, α such that item 3 of Assumption 1 holds. Use the algorithm presented in Section 2.1 of [13], obtaining the boolean v and the set of realsolutions R . If R is the empty set, set u as true, else set u as false. Output the booleans v and u , and the set R . Remark 3.
An appropriate choice of y ∈ R N \ V R ( (cid:98) g ε ) could be p itself, or p + N (0 , σ ) for some randommultivariate Gaussian noise with mean zero and variance σ . In step 4 above, if items 2, 3, 4, or 5 ofAssumption 1 fail to hold, then generating new and random points z, γ, y, α could be required.8 xample 1. As an illustrative example, we would like a configuration which we know to be locally rigid,but fails to be infinitesimally rigid so that Theorem 1 does not apply, and which is also a singular pointof the member constraints so that the local real dimension may differ from the local complex dimensioncomputed in [30]. As a simple example, consider a graph on nodes [5] = { , , , , } with edges E = { , , , , , , } . Consider the configuration in Figure 4. As can be verified by calculating theFigure 4:nullspace of the Jacobian at generic points and at this point, this configuration is singular. Its Jacobianadmits an additional infinitesimal flex, and hence Theorem 1 does not apply. However, it is also simpleenough that we can see it is locally rigid: Since the triangles among nodes 123 and 124 are both rigid, theonly node that could possibly move in a continuous flex is node 5. However, we can also see that node 5 isrestricted by node 3 to move in one circle of radius equal to the edge length of edge 35, while by node 4 itis restricted to another similar circle. These circles intersect in exactly one point, the location of node 5 inour initial configuration p . Thus p is locally rigid.We implemented Algorithm 1 for this example using HomotopyContinuation.jl in julia , obtaining ε -local rigidity for ε ∈ { . , . , . , . } . The point y was generated by Gaussian noise applied tothe original configuration p , then paths were tracked, obtaining zero real solutions. For each value of ε all paths were trackable and the items in Assumption 1 were satisfied, but none of the resulting solutionswere real-valued. Therefore, we can conclude by Theorem 5 that this configuration p is ε -locally rigid for ε ∈ { . , . , . , . } . Other choices of ε can be made by another user according to their application.Figure 5: Several singular configurationsThe example we just described is a singular configuration depicted in Figure 5 (left). Now consider thesingular configuration depicted in Figure 5 (right). This configuration is clearly not locally rigid. However,again it is singular and admits infinitesimal flexes, hence Theorem 1 does not apply. In this singular config-9ration, nodes 3 and 4 coincide, yielding a flexible node 5 which can freely move around a circle centered atnodes 3 ,
4. This configuration is obviously flexible, and when we run Algorithm 1 we obtain exactly two realsolutions (after projecting away the λ components). Upon examination, these two real solutions correspondto node 5 moving in either possible direction along a circle centered at nodes 3 ,
4, as expected. Algorithm 2described in Section 5 may be applied to visualize this flex.
Example 2.
Consider the 3-prism of Figure 1. We have nodes [6] = { , , , , , } and edges E = { , , , , , , , , , , , } . As mentioned previously, for generic configurations the 3-prism isinfinitesimally rigid and therefore rigid. However, for the configuration p (cid:55)→ (cid:98) p p = p p p p p p p p p p p p p p p p p p = − √ − − √ − √ − √ −
30 1 3 (cid:55)→ . . . . . . . − . . . − . . − . − . . . . . (5)there is an infinitesimal flex (right side of Figure 1) which is linearly independent of the (cid:0) d +12 (cid:1) = 6 infinitesimalrigid motions. Therefore Theorem 1 fails to apply.Applying Algorithm 3 to gather some preliminary numerical evidence for how this configuration mightdeform, we obtain many approximate, numerical configurations near p , including, for example . . . . . . . − . . . − . . − . − . . . . . . (6)This suggests the configuration p can twist downwards, since the z coordinate of nodes 4 , , . . p which exactly satisfies the member constraints.Therefore we apply Algorithm 1 to p . A naive total degree homotopy [27, Section 8.4.1] would requiretracking 67 , ,
864 paths. However, using mixed volumes of Newton polytopes [20] lowers the path-trackingcount to just 1 , , julia coderequired is available at [16], which also includes code for an expository article [17] which treats the 3-prismexample in detail. In this calculation, we obtain ε -local rigidity for ε = 0 .
1, even though the computedconfigurations like (6) are outside such a small ε -ball. Therefore, there can be no real-valued flex connecting p and (6). In fact, since the 3-prism is famous, the rigidity of this initial configuration p is well-known,thus our calculations agree with previous results. In this section we describe an algorithm that repeatedly solves a system of parametrized polynomial equationsin order to produce a sequence of real-valued, valid configurations. If a continuous flex of p exists, then thisprocedure will produce a discrete sampling of points from that continuous flex. The resulting sequence ofconfigurations may be plotted and animated, yielding easily understandable information for the scientist. In10uture work, we plan to implement this algorithm in a freely available julia package, utilizing the existingalgorithms of the julia package HomotopyContinuation.jl [7], which implements polynomial homotopycontinuation as discussed in Section 2. A main goal for our package will be ease of use. This is currentlyunder development, but a rough example of its output is shown in Figure 6. These are images of a discreteFigure 6: 100 configurations of a cube flexingflex with M = 100 points computed using homotopy continuation on a cube deforming freely. The cube isobviously flexible, and the goal is to implement software which will find more surprising flexes from other,more complicated examples. We also note that other excellent software exists for homotopy continuation,including [4, 22, 29]. Algorithm 2:
A discrete flex
Input:
Initial configuration p ∈ R nd , edge set E , choice of ε > M ∈ Z > . Result:
A discrete flex in the form of a list P of configurations p , p , . . . , p M to be animated andvisualized, or potentially the message of ε -local rigidity for some ε = j · ε . Initialize a list with one element P = [ p ] to be filled with more points p j ∈ R N as in p , p , . . . , p M ifthe algorithm succeeds. for j in 1:M do Set ε := j · ε . Apply Algorithm 1 with inputs p , E, ε , collecting the outputs u, v, R . if u = 1 then Output the current list P and the message that p is ε -locally rigid. else Collect the output set R and store it as R j = R . Set p j := argmin { dist( p j − , w ) : w ∈ R j } . Append p j to the list P . end end Return the list P as well as an animation of each of its M configurations displayed in R d if d = 2 , p . Consider the square system of11quations F y,ε ( x , . . . , x N , λ , λ ) = (cid:98) g ε ( x ) λ ( x − y ) + λ ∇ (cid:98) g ε α λ + α λ − : C N +2 → C N +2 . After solving and finding a new configuration p at radius ε away from p , this means we have obtainedsolutions ( x, λ ) to F y,ε for some specific y ∈ R N and ε >
0, where one of these solutions has x = p . Wecould then consider the homotopies perturbing y to y (cid:48) or perturbing ε to ε (cid:48) as in H y ( x, λ, t ) = F (1 − t ) y (cid:48) + ty,ε ( x, λ ) (7)or H ε ( x, λ, t ) = F y, (1 − t ) ε (cid:48) + tε ( x, λ ) , (8)either of which would produce new and relevant configurations. In particular, using H ε can allow us togenerate new solutions for expanding (or contracting) ε -balls about p , slightly modifying line 4 of Algorithm2. Remark 4.
In the homotopies above we have removed the usual factor γ . This allows real-valued solutionslike p to stay real-valued along the parameter homotopy. This will succeed unless we cross the discriminant,so we can expect success for | ε − ε (cid:48) | small. For small perturbations ε to ε (cid:48) we can numerically follow acontinuous flex and every step of the path-tracking process will compute new and real-valued configurationssampled from that continuous flex. More precisely, starting from a non-singular solution to F y,ε we cancontinue to a non-singular solution of F y,ε (cid:48) for | ε − ε (cid:48) | small. Even if we do cross the discriminant, it maybe due to phenomena involving complex solutions elsewhere, and thus not affect our particular real-valuedsolution. Remark 5.
In particular, letting ε → p on the ε -spheretowards p . If there is a continuous flex of p we can expect one of the real-valued solutions on the ε -sphereto move towards p as ε →
0. The collection of points computed along the way are a discrete flex of p ,having been sampled from a continuous flex. We can also follow the discrete flex away from p by letting ε increase, computing points in a parameter homotopy as above. Remark 6.
It is tempting to use monodromy or the trace test on each of the p j to ensure they are points onthe same irreducible component. However, even if they are on the same irreducible component, they could beon distinct connected components of the real algebraic set. Therefore, applying the trace test or monodromyare necessary but not sufficient conditions for our discrete flex to be sampled from the same connectedcomponent of V R ( (cid:98) g ). This option could be added to Algorithm 2 explicitly, but this would require witnesssets for each irreducible component of V ( (cid:98) g ) be computed, whereas here we are only assuming computationof a witness set for the pure ( N − V ( (cid:98) g ε ).Here we will describe another useful application of homotopy continuation to the problem of rigidity.Presented with a new example, this would be the first algorithm to try. The output of Algorithm 3 givesinformation about how the example may deform, should flexes exist, and can be applied even to much largerexamples. The only requirement is that you calculate an infinitesimal flex, since Algorithm 3 specificallysearches in that direction of configuration space. 12 lgorithm 3: Infinitesimal flex parameter homotopy
Input:
Initial configuration p ∈ R nd , edge set E , infinitesimal flex v ∈ R nd , choice of ε > M ∈ Z > . Result:
A discrete flex in the form of a list P of configurations p , p , . . . , p M to be animated andvisualized, or potentially the message of no real solutions for some ε = j · ε . Change coordinates to the moving frame as in Theorem 2. Let (cid:98) p and (cid:98) v be the resultingconfiguration and infinitesimal flex. Initialize a list with one element P = [ (cid:98) p ] to be filled with morepoints p j ∈ R N as in [ (cid:98) p , p , p , . . . , p M ] if the algorithm succeeds. for j in 1:M do Set ε := j · ε , set (cid:96) v = (cid:98) v T x − (cid:98) v T p , and set (cid:96) v,ε = (cid:98) v T x − (cid:98) v T p − ε . Set g as the polynomial system of moving frame member constraints associated to p as inDefinition 3. Solve the real parameter homotopy, collecting the resulting solution as p j by following thesolution p j − from t = 1 to t = 0 in h t ( x ) = (1 − t ) (cid:20) g(cid:96) v,(cid:15) (cid:21) + t (cid:20) g(cid:96) v (cid:21) . To be clear, p j − solves h ( x ). Follow that one solution via homotopy toward a new solutionstored as p j , which satisfies h ( x ). if p j is complex-valued or the parameter homotopy failed then Output the current list P and the message that there were no real solutions past ε = j · ε . else Collect the real-valued output p j and append it to the list P . end end Return the list P as well as an animation of each of its M configurations displayed in R d if d = 2 , Remark 7.
Here we discuss how to interpret the output of Algorithm 3. Although the configurationsobtained in Algorithm 3 are valid approximate solutions to the polynomial system up to standard numericaltolerances of your software, deciding whether a given point is real-valued or complex-valued is still typicallydone using a cutoff threshold. Therefore, the numerically computed configurations p , p , . . . , p M may trulycorrespond to complex-valued solutions with extremely small imaginary parts.However, they are not without meaning. Consider Examples 1, 2, and 3. For Example 1, Algorithm 3 didnot return any solutions even when we used the infinitesimal flex directions to search. This means the resultscorrectly suggest the configuration is very rigid. For both the 3-prism and Cluster 891529 in Examples 2 and3, whenever we used a randomly chosen vector w that was not pointing in the direction of any infinitesimalflex v of p , then Algorithm 3 correctly returned no output. Only when we specifically searched toward aninfinitesimal flex direction did the algorithm return real-valued solutions. Therefore, Algorithm 3 gives validand useful information about near-deformations , even if those numerically computed deformations turn outto be complex-valued configurations satisfying the member constraints but with extremely small imaginaryparts. This is still useful and relevant information.If further information is required of these numerically computed configurations p , p , . . . , p M , it is possibleto use Smale’s α -theory to determine if the p i are so-called approximate zeros of the given polynomial system[6, 26]. Roughly, this decides whether the computed configuration is within the region of quadratic conver-gence of Newton’s method to some exact solution of the system of polynomials. In fact, it is also possibleto decide if a given numerical solution corresponds to an exact real-valued solution. This is implemented inthe software alphaCertified [15]. Example 3.
We examine Cluster 891529 which is a so-called sticky sphere cluster. Briefly, nano andmicroscale particles sometimes interact only on distance ranges much smaller than their diameters, so itmakes sense to model them taking this into account. For more details see [18]. In this realm, Cluster 89152913as shared with us by Miranda Holmes-Cerfon via personal communication. It is believed that Cluster891529 is rigid but not second-order rigid (for the definition of second-order rigid see [9]). In particular, itis not first-order rigid (Theorem 1) since it admits infinitesimal flexes. We calculated an infinitesimal flex v Figure 7: Algorithm 3 applied to sticky sphere Cluster 891529numerically using the Jacobian of the system of member constraints, and then applied Algorithm 3 to thisinitial configuration p of Cluster 891529. The output is displayed in Figure 7. We find experimentally thatmany nearby numerical configurations can be obtained using v , or in fact any vector that points in roughlythe same direction in configuration space. This suggests that if Cluster 891529 is locally rigid, it is a singularpoint of of the system of member constraints of high multiplicity.Unfortunately, larger examples like Cluster 891529 of Example 3 are currently out of reach of Algorithm1 to test for ε -local rigidity. However, there is good news. HomotopyContinuation.jl version will soonbe finished, which will include significant changes in how polynomials are handled. In particular, insteadof expanding every polynomial in the monomial basis, the inherent structure present in polynomial systemsarising in applications will be preserved throughout, improving memory requirements, and the numerics offunction evaluations. For example, the polynomial systems we deal with here involve many sums of squares ofdifferences of variables. Expanding these in the monomial basis results in much higher memory requirements.When we square the entire polynomial itself, as in Definition 5, it is in fact much worse. There is good reasonto believe this example will be within reach once this new software is finished. Nonetheless, Cluster 891529clearly demonstrates the utility of Algorithm 3.
We considered bar frameworks from the perspective of numerical algebraic geometry, in particular usingresults from real algebraic geometry and polynomial homotopy continuation. We prove Theorem 5 whichsupports Algorithm 1 for testing the ε -local rigidity of a framework. To show a framework in configuration p ∈ R nd is ε -locally rigid is to show that any continuous deformations of p remain nearby to p inconfiguration space. In particular, they do not exit a sphere of radius ε about p in configuration space.14n ε -locally rigid framework may be locally rigid or flexible (see Figure 2), but since ε can be freely chosenby the user, this makes ε -local rigidity a useful property in any application where very small movementscan be ignored. Another advantage is that ε -local rigidity applies to any configuration p , be it smooth orsingular, generic or non-generic. Finally, we present Algorithms 2 and 3 which also use polynomial homotopycontinuation to provide useful visual information for the scientist, even in examples that are currently toolarge for Algorithm 1. References [1] Timothy Abbott, Reid Barton, and Eric Demaine. Generalizations of Kempe’s Universality Theorem.
Master’s Thesis, MIT , 2009.[2] L. Asimow and B. Roth. The rigidity of graphs. II.
J. Math. Anal. Appl. , 68(1):171–190, 1979.[3] Philippe Aubry, Fabrice Rouillier, and Mohab Safey El Din. Real solving for positive dimensionalsystems.
J. Symbolic Comput. , 34(6):543–560, 2002.[4] Daniel J. Bates, Jonathan D. Hauenstein, Andrew J. Sommese, and Charles W. Wampler. Bertini:Software for numerical algebraic geometry. Available at bertini.nd.edu with permanent doi:dx.doi.org/10.7274/R0H41PB5.[5] Daniel J. Bates, Andrew J. Sommese, Jonathan D. Hauenstein, and Charles W. Wampler.
NumericallySolving Polynomial Systems with Bertini . Society for Industrial and Applied Mathematics, Philadelphia,PA, 2013.[6] Lenore Blum, Felipe Cucker, Michael Shub, and Steve Smale.
Complexity and real computation .Springer-Verlag, New York, 1998. With a foreword by Richard M. Karp.[7] Paul Breiding and Sascha Timme. Homotopycontinuation.jl: A package for homotopy continuation injulia. In James H. Davenport, Manuel Kauers, George Labahn, and Josef Urban, editors,
MathematicalSoftware – ICMS 2018 , pages 458–465, Cham, 2018. Springer International Publishing.[8] R. Connelly and S.D. Guest.
Frameworks, Tensegrities and Symmetry . To appear, 2020.[9] Robert Connelly. The rigidity of certain cabled frameworks and the second-order rigidity of arbitrarilytriangulated convex surfaces.
Adv. in Math. , 37(3):272–299, 1980.[10] Robert Connelly and Herman Servatius. Higher-order rigidity—what is the proper definition?
DiscreteComput. Geom. , 11(2):193–200, 1994.[11] Grigore Gogu. Mobility of mechanisms: a critical review.
Mech. Mach. Theory , 40(9):1068–1097, 2005.[12] Jacob E. Goodman, Joseph O’Rourke, and Csaba D. T´oth, editors.
Handbook of discrete and compu-tational geometry . Discrete Mathematics and its Applications (Boca Raton). CRC Press, Boca Raton,FL, 2018. Third edition of [ MR1730156].[13] Jonathan D. Hauenstein. Numerically computing real points on algebraic sets.
Acta ApplicandaeMathematicae , 125(1):105–119, Sep 2012.[14] Jonathan D. Hauenstein and Andrew J. Sommese. What is numerical algebraic geometry?
Journal ofSymbolic Computation , 79:499 – 507, 2017. SI: Numerical Algebraic Geometry.[15] Jonathan D. Hauenstein and Frank Sottile. Algorithm 921: alphaCertified: certifying solutions topolynomial systems.
ACM Trans. Math. Software , 38(4):Art. 28, 20, 2012.[16] Alexander Heaton. Tensegrity. https://github.com/alexheaton2/tensegrity . Accessed: 2020-03-01. 1517] Alexander Heaton. Nonlinear algebra via tensegrity structures. arXiv e-prints , page arXiv:1908.08392,August 2019.[18] Miranda Holmes-Cerfon. Sticky-sphere clusters.
Annual Review of Condensed Matter Physics , 8(1):77–98, 2017.[19] Miranda Holmes-Cerfon, Louis Theran, and Steven J. Gortler. Almost-rigidity of frameworks, 2019.[20] Birkett Huber and Bernd Sturmfels. A polyhedral method for solving sparse polynomial systems.
Math.Comp. , 64(212):1541–1555, 1995.[21] Sergi Hern`andez Juan and Josep M. [Mirats Tur]. Tensegrity frameworks: Static analysis review.
Mechanism and Machine Theory , 43(7):859 – 881, 2008.[22] T. L. Lee, T. Y. Li, and C. H. Tsai. HOM4PS-2.0: a software package for solving polynomial systemsby the polyhedral homotopy continuation method.
Computing , 83(2-3):109–133, 2008.[23] Peter Olver. Lectures on moving frames. , 2012.[24] F. Rouillier, M.-F. Roy, and M. Safey El Din. Finding at least one point in each connected componentof a real algebraic set defined by a single equation.
J. Complexity , 16(4):716–750, 2000.[25] A. Seidenberg. A new decision method for elementary algebra.
Ann. of Math. (2) , 60:365–374, 1954.[26] Steve Smale. Newton’s method estimates from data at one point. In
The merging of disciplines: newdirections in pure, applied, and computational mathematics , pages 185–196. Springer, New York, 1986.[27] Andrew J Sommese and Charles W Wampler.
The Numerical Solution of Systems of PolynomialsArising in Engineering and Science . World Scientific, 2005.[28] G. Strang.
Computational Science and Engineering . Wellesley-Cambridge Press, 2007.[29] Jan Verschelde. Algorithm 795: Phcpack: A general-purpose solver for polynomial systems by homotopycontinuation.
ACM Trans. Math. Softw. , 25(2):251–276, June 1999.[30] Charles W. Wampler, Jonathan D. Hauenstein, and Andrew J. Sommese. Mechanism mobility and alocal dimension test.
Mechanism and Machine Theory , 46(9):1193 – 1206, 2011.[31] Walter Whiteley. Matroids and rigid structures. In
Matroid applications , volume 40 of