Dorst L., Fontijne D., Mann S. Geometric Algebra for Computer Science. An Object Oriented Approach to Geometry

Подождите немного. Документ загружается.

FUNDAMENTAL

ALGORITHMS FOR

NONLINEAR PRODUCTS

In the previous chapter, we looked at how to implement the linear products in geometric

algebra. The linearity of these products allowed us to implement them using linear alge-

bra or through a simple double loop. However, there are other operations in geometric

algebra that are nonlinear (such as inverse,

meet, join, and factorization). These cannot

be implemented in the same way.

In this chapter, we discuss the implementation of such nonlinear geometric algebra oper-

ations. The nonlinearity results in more complex algorithms, still reasonably efﬁcient but

typically an order of magnitude more time-consuming than linear operations.

We give algorithms for the inverse, for exponentiation, for testing whether a multivector

is a blade or a versor, for blade factorization, and for the efﬁcient computation of

meet

and join.

21.1 INVERSE OF VERSORS (AND BLADES)

We need to compute the inverse of the elements we construct in geometric algebra. Those

are almost exclusively blades or versors (the only exceptions were the bivectors in an

exponent, and we have no need to invert those). Invertible blades are always versors, since

529

530 FUNDAMENTAL ALGORITHMS FOR NONLINEAR PRODUCTS CHAPTER 21

a k-blade (the outer product of k vectors) can always be written as the geometric product

of k vectors (i.e., a k-versor); see Section 6.4.1. So if we can invert a versor V, we can also

invert blades and thus most of the elements we construct.

There is a very efﬁcient way to invert versors, which we call the versor inverse method

(because, in general, it works only on versors). This inversion is almost trivial. It starts

with the observation that the quantity V



V is always scalar for any versor V. This is easily

seen by factorizing the versor:



V = (v

···v

)(v

···v

) = (v

···(v

) v

) ···v

) = v

···v

The inverse of the versor is then simply

−1



V/(V



V ).

(21.1)

This method would fail if makes V



V = 0. But in that case (which may occur when one

of the vector factors is a null vector), the multivector V is actually not a versor (since all

versors are invertible by deﬁnition). V is then not actually a versor since versors are always

invertible.

21.2 INVERSE OF MULTIVECTORS

As long as you stay within the bounds of geometric algebra (as opposed to Clifford algebra;

see Section 7.7.2 for our take on the difference), the versor inverse method should cover

all of your inversion needs. But should the need arise to invert a general multivector, you

can use the following technique to invert any invertible multivector. It is based on the

linear algebra implementation approach described in Section 20.1. Since

 [[ A

]] [[ B]] ,

it follows that

−1

B  [[ A

]]

−1

[[ B]]

(just left-multiply both sides by A and [[ A

]] ). So the matrix of the inverse of A is the

inverse of the matrix of A. This insight can be used to invert any multivector A as long as

the matrix [[ A

]] is invertible, using the following steps:

•

Compute the geometric product matrix [[ A

]] ;

•

Invert this matrix;

•

Extract [[ A

−1

]] as the ﬁrst column of [[ A

]]

−1

To understand why the last step is correct, inspect the algorithm for computing the geo-

metric product matrix. With the scalar 1 as the ﬁrst element of

L, the ﬁrst column of the

matrix [[ A

]] is constructed from the geometric product of A with the ﬁrst element of L,

which is 1. Therefore the ﬁrst column of a geometric product matrix [[ A

]] is [[ A]] itself

(this may also be observed in Figure 20.1). Hence we ﬁnd the coefﬁcients of A

−1

as the

ﬁrst column of [[ ( A

−1

)

]] = [[ A

]]

−1

SECTION 21.3 EXPONENTIAL, SINE, AND COSINE OF MULTIVECTORS 531

This inversion method has two main disadvantages: it is both slow and numerically

imprecise due to ﬂoating point round-off errors. Therefore, you should not use it to invert

a blade or a versor—use the method of Section 21.1 for those.

21.3 EXPONENTIAL, SINE, AND COSINE

OF MULTIVECTORS

The exponential and trigonometric functions were introduced for blades in Section 7.4.2.

For easy reference, we repeat the equations for the polynomial expansion of the

exponential:

exp(A) = 1 +

+ ···

If the square of A is a scalar, this series can be easily computed using standard trigono-

metric or hyperbolic functions:

exp(A) =

⎧

⎪

⎨

⎪

⎩

cos α + A

sin α

if A

= −α

1 + A if A

= 0

cosh α + A

sinh α

if A

= α

where α is a real scalar. Since both a versor and a blade have scalar squares, their exponen-

tials are easy to compute by directly applying these formulas.

Unfortunately, we will need to take exponentials of nonblade elements. The exponen-

tials of general bivectors generate the continuous motions in a geometry, and only in

fewer than four dimensions are bivectors always 2-blades. In the important conformal

model, rigid body motions are exponentials of bivectors that are not 2-blades, as we saw

in Section 13.5. So we need techniques to process them.

When the special cases do not apply, the series for the exponential can be evaluated explic-

itly up to a certain order. This tends to be slower and less precise. Experience shows that

evaluating the polynomial series up to order 10 to 12 gives the best results for 64 bit dou-

bles. For the exponential, a rescaling technique is possible that will increase accuracy, as

follows.

Suppose you want to compute exp(A).IfA is a lot larger than unity, A

can be so large that

the series overﬂows the accuracy of the ﬂoating-point representation before it converges.

To prevent this problem, we can scale A to near unity before evaluating the series, because

of the identity

exp(A) =



exp







For our purposes, any s ≈A will do. In a pra ctical implementation, we choose s to be a

power of two, so that we can efﬁciently compute



exp



by repeatedly squaring exp(

532 FUNDAMENTAL ALGORITHMS FOR NONLINEAR PRODUCTS CHAPTER 21

21.4 LOGARITHM OF VERSORS

We do not know of a general algorithm for computing the logarithm of arbitrary versors.

However, for many useful cases a closed-form solution was found; see, for example,

Section 10.3.3 (rotation), Section 13.5.3 (rigid body motion), and Section 16.3.4 (pos-

itively scaled rigid body motion). We have not yet included these logarithms in our refer-

ence implementation, but you can ﬁnd a C++ implementation of these logarithms in the

GA sandbox source code package.

21.5 MULTIVECTOR CLASSIFICATION

At times, it may be useful to have an algorithm that classiﬁes a multivector as either a

blade, a versor, or a nonversor. For instance, classiﬁcation is useful as a sanity check of

interactive input, or to verify whether and how a result could be displayed geometrically.

Yet testing whether a multivector is a versor or a blade is nontrivial in our additive

representation. A blade test requiring that the multivector is of a single grade is insufﬁcient

(for example, e

∧ e

+ e

∧ e

is of unifor m grade 2, but not a blade). Adding the rule that

the square of the multivector must be a scalar does not help (since e

∧e

+ e

∧e

squares to −2 in a Euclidean metric, but it is not a blade).

The classiﬁcation algorithm below (from [7]) performs this test correctly. It is one algo-

rithm that can be used either for the versor test or for the blade test. It is natural that these

tests should structurally be very similar, for all invertible blades are also versors. As you

glance through it, you notice that the algorithm uses inverses. This is ﬁne for a versor test,

since versors need to be invertible by deﬁnition. However, a blade need not be invertible,

and yet null blades are still blades.

In fact, “being a blade” is deﬁned as “factorizable by the outer product”, and that

does not depend on a metric at all. Our reference [7] makes good use of this freedom:

if you are testing whether V is a blade, you must perform the geometric products in

the algorithm using a Euclidean metric. This choice of convenience eliminates all null

blades that would have had to be taken into account without actually affecting factor-

izability, and that simpliﬁes the algorithm. (As an example, the multivector e

∧∞ is

a blade, but in the conformal metric contains a null factor ∞. This makes the blade

noninvertible, but it is of course still a blade. By using a Euclidean metric, ∞ is treated

as a regular vector, the blade becomes invertible, and the test can run as it would for

∧ e

On the other hand, determining whether a multivector V is a versor (i.e., a geometric

product of invertible vectors) clearly depends on the precise properties of the metric.

So for a versor test, you have to run this algorithm in the actual metric of the algebra

of the versor.

SECTION 21.6 BLADE FACTORIZATION 533

The classiﬁcation process for V consists of three parts:

1. Test if the versor inverse



V/(V



V ) is truly the inverse of the multivector V. This

involves the following tests:

grade(



−1

)

= 0,



−1

= V

−1



If either of these test fails, we can report that the multivector is a nonversor, and

hence a nonblade.

The use of grade involution in the two equations above is an effective trick from [9].

It prevents multivectors that have both odd and even grade parts from sneaking

through the test. If V is an even versor, then



V = V.IfV is an odd versor,



V = −V.

But if V has both odd and even grade parts, the grade involution prevents odd and

even parts from recombining in such a way that they cancel each other out. See

structural exercise 21.8 for an illuminating example.

2. The second test is on the grade preservation properties of the versor: applying a

versor (and hence an invertible blade) to a vector should not change the grade. So

for each basis vector e

of the vector space R

, the following should hold if V is a

versor:

grade(



V e



V )

= 1.

When the multivector does not pass this test, we report a nonversor (and hence a

nonblade); otherwise we know that it is either a versor or a blade.

3. The ﬁnal part makes the distinction between blades and versors by simply checking

whether the multivector is of a single grade. If so, it is a blade; otherwise it is a

versor.

Unfortunately, we currently have no general solution to the problem of determining

the ty pe of multivectors that are not quite blades or not quite versors in a numerically

informed manner. We would like to correct any numerical drift in these fundamental

properties due to repeated ﬂoating point round-off errors, but we do not know how to

make versors or blades out of almost-versors and almost-blades in an optimal way in gen-

eral metrics.

21.6 BLADE FACTORIZATION

This section deals with generic factorization of blades. That is the problem to ﬁnd, for a

given blade B of grade k,asetofk vectors b

such that

B = b

∧ b

∧···b

You may want to factorize a blade because you want to use it as input to libraries that

cannot handle blades (such as standard libraries for linear algebra or computer graphics),

534 FUNDAMENTAL ALGORITHMS FOR NONLINEAR PRODUCTS CHAPTER 21

so that you are required to process their vector factors separately. Or you may need a

factorization to implement another low-level algorithm such as the computation of

meet

and join below.

We are only concerned with the outer product and consequently (as in the previous

section) are allowed to choose any convenient metric. To avoid any problems with null

vectors, we use the Euclidean metric for all metric products in the algorithm that follows.

A ﬁrst step towards a useful factorization algorithm is ﬁnding potential factors of a blade

B. This can be done using projection. Take any candidate vector c (for example, a basis

vector), and project it onto B:

f = (cB) B

−1

If f is 0, then another candidate vector should be tried. For better numerical stability, we

should also try another candidate when f is close to 0 (or better yet, try many different

candidate vectors and use the one that results in the largest f ). In any case, when f is a

usable factor we can remove it from B simply by dividing it out:

= f

−1

B.

(21.2)

Then we can write (see structural exercise 2)

B = f ∧ B

so that we have found our ﬁrst factor of B. We now repeat this process on the blade B

iteratively, until we are left with a ﬁnal vector factor.

This basic idea works, but the procedure may be inefﬁcient in high-dimensional spaces.

Many candidate vectors may result in a zero value for f, which is a waste of effort. It is

better to use the structure of B itself to limit the search for its factors in an efﬁcient, stable

technique, as follows. The blade B, as given, is represented as coordinates relative to a basis

of blades. We take the basis blade with the absolute largest coordinate on this basis. Let

that be E = e

∧ e

∧ ... ∧ e

. We then use the basis vectors that make up E as candidate

vectors for projection. This selection procedure guarantees that the projection of each of

the candidate vectors is nonzero (we show this at the end of this section).

One ﬁnal issue is the scaling of each of the factors. Due to the projection, the factors

do not have a predictable scale, which is an awkward property. Our implementation of

the factorization algorithm normalizes the factors and returns the scale of the blade as a

separate scalar.

1 We remove the factor f from B, because we want the factors to be orthogonal. If you just want factors that make

B and don’t care about orthogonality or scale, you may adjust the algorithm to use the original B in each loop

of the ﬁnal algorithm below. But you must then make sure that your factors are linearly independent.

2 Other solutions are to premultiply the ﬁrst factor with the scale, or to apportion the scale evenly over each factor.

Which method is most convenient may partly depend on the subsequent use of the factorization. In any case,

the unit factors with a separate scale can be transformed into any of the other representations easily.

SECTION 21.6 BLADE FACTORIZATION 535

The ﬁnal algorithm [6, 8] for factorizing a blade becomes:

1. Input: a nonzero blade B of grade k.

2. Determine the norm of B: s = B.

3. Find the basis blade E in the representation of B with the largest coordinate;

determine the k basis vectors e

that span E.

4. Let the current input blade be B

← B/s.

5. For all but one of the basis vectors e

of E:

(a) Project e

onto B

: f

= (e

B

) B

−1

(b) Normalize f

. Add it to the list of factors.

: B

← f

−1

B

6. Obtain the last factor: f

= B

. Normalize it.

7. Output: the factors f

and the scale s.

Some notes on metric are in order. First, the algorithm also works for null blades, since no

blade is actually null in the Euclidean metric that is used during the factorization algo-

rithm itself. Second, the output of the algorithm is a set of orthonormal factors in the

Euclidean metric that was used within the algorithm. That may not be the metric of the

space of interest. If you desire orthonormality in some other metric, construct the metric

matrix of the factors, perform an eigenvalue decomposition, and use this to construct an

orthonormal set of factors (see also Section 19.4).

To see why the projection of e

on B

is never zero (step 5a), note that there always exists

arotorR that turns E to the original B. This rotor will never be over 90 degrees (for that

would imply that EB = 0, yet we know that EB must be nonzero to be the basis blade

with the largest coordinate in B). We may not be able to compute the rotation R easily, but

we can ﬁnd R

= RR= BE

−1

. Since R is never over 90 degrees, R

willneverbeover180

degrees. Because of this, the quantity

+ R

) must be nonzero: no e

is rotated

far enough by R

to become its own opposite. We rewrite this and ﬁnd:

0 =

+ R

)

+ BE

−1

)

− BE

−1



−1

)

−



−1

)

= (e

B) B

−1

Therefore none of the e

from E projects to zero on B,sof

is nonzero in the ﬁrst pass of

the algorithm. After removal of f

, the same argument can be applied to the next e

on B

and so on. Hence none of the f

are zero.

536 FUNDAMENTAL ALGORITHMS FOR NONLINEAR PRODUCTS CHAPTER 21

21.7 THE MEET AND JOIN OF BLADES

When we introduced the meet and join of blades in Chapter 5, we saw that they are in

some sense the geometrical versions of intersection and union from set theory. We can

use this correspondence to compute them. To simplify the description of the algorithm,

we ﬁrst illustrate some of the basic ideas with regular set theory and Venn diagrams and

then transfer them to

meet and join.

Figure 21.1(a) shows a Venn diagram of two nondisjoint sets A and B, their union, A ∪ B,

and their intersection, A ∩ B. We introduce a symmetric set difference through a delta

product, deﬁning A

M B as A ∪ B minus A ∩ B. T his is illustrated in Figure 21.1(d). We

also use the complement (A

M B)

∗

as illustrated in Figure 21.1(e). The dashed line along

the border of Figure 21.1(e) indicates that ( A

M B)

∗

extends to include all elements that

are not in A or B.

Suppose that we have a function s() that determines the size of a set. We can relate the size

of the various sets through

s(A ∪ B) =

s(A) + s( B) + s(A

M B)

. (21.3)

To see why this holds, ﬁrst convince yourself that

s(A) + s( B) = s(A ∪ B) + s(A ∩ B),

and then superimpose Figure 21.1(b), Figure 21.1(c), and Figure 21.1(d). Note that every

area is covered twice, hence the division by 2 in (21.3). Likewise,

A ∪ B

A ∩ B

(A  B)*

A  B

Figure21.1: Venn diagrams illustrating union, intersection, and the delta product of two sets.

(a) Two sets A and B, (b) the union A ∪ B, (c) the intersection A ∩ B, (d) the delta product A

 B,

and (e) The dual of the delta product (A

 B)

∗

SECTION 21.7 THE MEET AND JOIN OF BLADES 537

s(A ∩ B) =

s(A) + s( B) − s(A

M B)

. (21.4)

The

meet and join products are the geometrical versions of set intersection and union,

respectively, applied to blades rather than sets. They really work on the bases spanning the

blades, and are in that sense discrete. We want to compute the

join of two blades A and

B. Instead of thinking of A and B as continuous oval blobs, think of them of as discrete

factors, united by a boundary. That converts Figure 21.1 into Figure 21.2. Let us assume

that A and B can be factored as

A = a

B = cb

This means that A and B have one factor c in common, and a

, a

, and b

are independent

(how we arrive at this factorization is discussed later on). Figure 21.2(a) illustrates this

basic setting. Obviously, A ∩ B is proportional to c (Figure 21.2(c)), and the A ∪ B is

proportional to a

∧ a

∧ c ∧ b

(Figure 21.2(b)). Figure 21.2(d) illustrates A M B, and

Figure 21.2(e) illustrates (A

M B)

∗

The algorithm for computing the

meet and join [6, 8] presented in this section ﬁrst

computes the required grade of the

meet and join, for which we need the delta product

for blades (from [8]). Then it works towards constructing either the

meet or the join.

As we observed in Section 19.3, the geometric product acts like a spatial exclusive or

on basis blades. The same is true for regular blades when we focus on highest nonzero

grade part of the geometric product: it is the spatial exclusive or of the blades since

it contains all the factors that do not occur in both arguments. In our example,

A ∪ B

A  B

(A  B)*

A ∩ B

Figure 21.2: Venn diagrams illustrating Meet, Join, and the delta product of two blades.

A = a

∧ a

∧ c, B = c ∧ b

(see also Figure 21.1).

538 FUNDAMENTAL ALGORITHMS FOR NONLINEAR PRODUCTS CHAPTER 21

AB= (a

c)(cb

) = (c · c) a

Note that we do not need to factorize the blades explicitly to get this result. The geometric

product automatically eliminates the dependent factors. So the top grade part of a geo-

metric product ABis a blade that contains all factors that were present in either A or B.

We call this the delta product of blades, and deﬁne it as follows:

M B ≡AB

max

where max is the largest grade such that A M B is not zero. The delta refers to the

property of the delta product to compute the blade that contains factors that are

present in either of the arguments, which makes it like the symmetric difference of

the factors in ABwhich we discussed above.

Once we have the delta product of blades, we can derive the following useful relations by

analogy to (21.3) and (21.4):

grade(A ∪ B) =

grade(A) + grade(B) + grade(A

M B)

(21.5)

grade(A ∩ B) =

grade(A) + grade(B) − grade(A

M B)

(21.6)

Our

meet algorithm starts with a scalar, and expands it by the outer product with new

vectors until it arrives at the true

meet. Potential factors of the meet satisfy the following

conditions, which you may verify using Figure 21.2.

•

They are not factors of the delta product;

•

They are factors of A and B.

Likewise, the algorithm initially assumes that the

join is the pseudoscalar, and removes

factors from it until the true

join is obtained. Factors that should not be in the join sat-

isfy the following conditions:

•

They are factors of the dual of the delta product;

•

They are not factors of A and B.

Before we can formulate the algorithm, we repeat two relations from Section 5.3:

A ∪ B = A ∧ ((A ∩ B)

−1

B)

(21.7)

A ∩ B = (B(A ∪ B)

−1

)A (21.8)

These equations can be used to obtain the

join from the meet, and vice versa. We now

have enough building blocks to construct an algorithm to compute the

meet and join:

1. Input: two blades A, B, and possibly a threshold  for the delta product.

2. If grade(A) > grade(B), swap A and B. This may engender an extra sign, so be

careful when you need to interpret the results (see Section 5.6).