Bayro-Corrochano E., Scheuermann G. Geometric Algebra Computing: in Engineering and Computer Science

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

458 D. Fontijne and L. Dorst

deﬁnes the join as the outer (wedge) product (i.e., join(A, B) =A ∧B) and the meet

as meet(A, B) = (A

∗

∧B

∗

)

∗

, where the dual

∗

is with respect to the whole vector

space. As a result, software packages based on these deﬁnitions (such as Bige-

bra/CLIFFORD [1]) require the dimension of the vector space to be set (dim_V

:= ...) before computing the meet. For example, in CLIFFORD:

> dim_V:=3: M3:=meet(e1we2,e2we3);

M3 := -e2

> dim_V:=4: M4:=meet(e1we2,e2we3);

M4 := 0

This shows that in CLIFFORD the meaning of intersection (

meet) depends on the

dimension of the full space, which we ﬁnd rather strange. So once again, we would

like to stress that this paper is about computing the

join (and indirectly the meet)of

blades in any geometric (possibly degenerate) situation. That is, the

join is the true

geometric subspace union, e.g.,

join(e

, e

) =e

and not join(e

, e

) =0.

In this paper, we present a new blade factorization algorithm for blades rep-

resented as a sum of based blades. We also use the new factorization algo-

rithm to construct a new algorithm to compute the

join. Both new algorithms

are improvements of previous work [2, 3, 6]. Our main contributions are the

FastFactorization algorithm, which factors blades by simply rearranging coordi-

nates, and the StableFastJoin algorithm, which efﬁciently computes the

join in a

numerically stable way.

Having a fast and stable

join implementation is especially important because it

encourages people to make proper use of geometric algebra. One of the great ad-

vantages of geometric algebra is that many equations are universal and that the

join

allows for universal subspace union. However, if it is perceived as slow in practice,

people will rewrite their equations to avoid the

join and replace it with linear inner

and outer products. This will often break the universality of their equations, splitting

one

join into many different cases (e.g., if the variable is a line, then do this; if the

variable is a plane, then do that; and so on).

We only brieﬂy consider computing the

meet (true subspace intersection) in

Sect. 5.3. Experimentation showed that adjusting our FastJoin algorithm to com-

pute the

meet directly is somewhat slower than computing it from the join using [6]

meet(A, B) =



Bjoin(A, B)

−1



A (1)

and also leads to more generated code. This is why the topic is ignored in the main

text.

We assume a Euclidean metric in all computations, as both factorization and

join

are metric-independent operations (i.e., we use the LIFT described in [4]). Through-

out the paper, n is the dimension of the vector space V

, and k is used to denote the

grade of the blade in the current context.

Efﬁcient Algorithms for Factorization and Join of Blades 459

2 Blade Factorization

The problem of factorizing a blade is the following. Given a k-blade B (in the addi-

tive representation), ﬁnd k vectors b

such that

B =b

∧b

∧···∧b

where for reasons of numerical stability, we prefer the factors b

to be “sufﬁciently

nonparallel.”

Not all k-vectors (homogeneous multivectors of grade k) are blades, i.e., not all

k-vectors can be factored (for example, e

∧ e

+ e

∧ e

). Such nonfactorizable

k-vectors are invalid input to our factorization algorithm. However, the algorithm

will not attempt to detect such input, since feeding it a nonblade does not have

catastrophic consequences. If required, a separate bladedness test as described in

[6] can be employed to make sure that the input to the algorithm is a blade.

To ﬁnd factors, one may project “probing vectors” p

onto the blade. In geometric

algebra, this is done by the projection operator [5, 6]



B

−1



B. (2)

If q

is not zero, it is a factor of B. By ﬁnding k linearly independent vectors q

a factorization of B has been found (up to scale). A straightforward choice for the

probing vectors p

are basis vectors e

. To obtain the factorization, one selects a total

of k independent projected vectors q

. This procedure is the essence of the outer

factorization algorithm presented in detail in [6], based on ideas in [3]. It works, but

the projection is computationally rather expensive.

2.1 New Algorithm for Blade Factorization

Figure 1 inspired our new factorization algorithm. The ﬁgure illustrates the steps

involved in the projection of a vector onto a 2-blade using geometric algebra by

means of the two contractions in (2). The ﬁrst inner product (p

B

−1

) computes the

orthogonal complement in B of the projected vector onto B (with a proportionality

factor of 1/B

). The second inner product in (p

B

−1

)B computes the orthog-

onal complement of the result within B, and this rotates the previous result in the

B-plane, to become the actual projection.

If we are only looking for factors, we are more interested in the fact that the ﬁnal

result q

is in B than in it being the precise projection of the original probing vec-

tor p

. However. this “being in B” is already guaranteed by the second contraction

in the projection equation. Hence, in our fast outer factorization algorithm, we save

on the ﬁrst contraction by using the orthogonal complement of the probing vector

with respect to the largest basis blade in B. The whole operation then becomes

merely the selection of appropriate coordinates. We turn this idea into an algorithm

as follows:

460 D. Fontijne and L. Dorst

Fig. 1 Orthogonal projection

of a vector onto a 2-blade,

displayed as the undoing of

the orthogonal-complement-

step of the inner product

Algorithm FastFactorization(B):

Let B be a k-blade with 1 <k<n(to exclude trivial cases). The algorithm computes

a factorization B =β b

∧b

∧···∧b

, where β is a scalar:

1. Find the basis blade F to which the absolute largest coordinate of B refers. Let f

be the basis vectors in F (i.e., F =f

∧f

∧···∧f

), and let β be the coordinate

that refers to F.

2. Compute B

=B/β.

3. For each f

, compute: b

=(f

F

−1

)B

The independency of the vectors b

guarantees that they form a factorization

of B

Theorem The factors b

as computed by the FastFactorization algorithm are lin-

early independent.

Proof B

is represented as a sum of grade k basis blades E

: B



(

)

j=1

. By

distributivity, each b



(

)

j=1

F

−1

)E

. Let us analyze the contribution of

each E

to each b

• If E

is equal to F, then (f

F

−1

)E

• Else if E

is orthogonal to F, we ﬁnd (f

F

−1

)E

=0.

• Else E

= σ

F

−1

) ∧ e

for some basis vector e

, where the sign σ

=±1

depends on the order of basis vectors in E

.Thise

cannot be equal to f

,for

then E

is equal to F, and also e

cannot be any of the other basis vectors in F,

for then E

would be 0 as it would contain the same basis vector twice. Thus, in

this case,



F

−1



E



F

−1







F

−1



∧e





F

−1



=±σ

where e

is not a factor of F.

Thus each b

equals



∧F=0

(±σ

)

Efﬁcient Algorithms for Factorization and Join of Blades 461

so that it does not contain any other factor of F than f

itself. Since the f

are linearly

independent, so are the b

. 

To illustrate the proof, let us compute the factors of the example B given above:

B =−0.8e

∧e

+0.4 e

∧e

−0.2 e

∧e

+0.6 e

∧e

Then F =e

∧e

and β =−0.8, and after scaling Bs becomes

∧e

−0.5 e

∧e

+0.25 e

∧e

−0.75 e

∧e

The factors are:



F

−1



B

= e

+0.25 e



F

−1



B

= e

+0.5 e

, (3)



F

−1



B

= e

−0.75 e

The diagonal typesetting of e

, e

should make it obvious that the b

are linearly

independent, and because the b

are all contained in B, and there are grade(B) of

them, they must form a factorization of B.

Even though the factors are linearly independent, they are not orthogonal in gen-

eral. For an estimation of how nonorthogonal the factors can be in the worst case,

let us assume that the B is placed so skewly relative to the basis that it has an equal

weight for each basis blade (it is not guaranteed that such an element B is indeed a

blade). For a unit blade B, this gives, as the worst case,

B =







(

)



j=1

±E

, so B

(

)



j=1

±E

For such an element B

, the largest possible absolute value of an inner product

between a pair of factors is the number of nonzero coordinates:

b

·b

≤n −k,

and the largest absolute value of an inner product between a pair of normalized

factors is



unit(b

) ·unit(b

)



≤

n −k

n −k +1

This last equation shows that the factorization method could be less usable in, for

example, 10-D space (e.g., if k =5 and n =10, then unit(b

) ·unit(b

)≤

) than

it is in, for example, 3-D space, because the factors are potentially less and less

orthogonal as the dimension of space n increases. However, in such spaces a factor-

ized blade representation becomes more attractive than the additive representation

462 D. Fontijne and L. Dorst

[10], and issues of efﬁcient factorization of blades in the additive representation are

then less crucial.

We still have to prove that the outer product of the computed factors b

actually

is equal to B

, i.e., that the scale and orientation of the computed factors are correct.

Theorem The outer product of the factors b

as computed by the FastFactorization

algorithm is equal to B

∧b

∧···∧b

Proof We know that the computed factors b

are actually in B

and are linearly

independent. So the proof can be simpliﬁed to proving that the largest (unit) basis

blade F of B

is correctly computed, because then the other basis blades will have

the right scale automatically . Now, F =f

∧f

∧···∧f

, and the contribution to

this basis blade of each computed factor b

is (f

F

−1

)F. Hence, if we can prove

that



F

−1



F =f

the full theorem holds. This is trivial because f

is contained in F, allowing us to

rewrite



F

−1



F =



−1



F =f



−1





3 Algorithms for Computing the Join of Blades

Given the FastFactorization algorithm of the previous section, it is straightforward

to formulate fast algorithms for the

join based on [2, 6].

3.1 Fast Join Algorithm

Given two blades A and B, the following algorithm computes the join J =join(A, B).

A small constant threshold value ε is required.

Algorithm FastJoin(A, B,ε):

1. Filter out trivial cases:

• If A and/or B is zero, return J =0.

• Else if A is of grade 0, return J =unit(B).

• Else if B is of grade 0, return J =unit(A).

• Else if A and/or B is of grade n, return J =I

• Otherwise continue with Step 2.

Efﬁcient Algorithms for Factorization and Join of Blades 463

2. If required, swap A and B so that grade(A) ≥ grade(B). If the swap is applied,

remember the “swapping-sign” σ = (−1)

grade(A) grade(B)

, otherwise set σ = 1.

The swap is for reasons of efﬁciency.

3. Set J ←unit(A).

4. Find the largest basis blade term F in B.

5. While grade(J) =n and not all basis vectors f

in F have been tried:

a. Take any basis vector f

in F which has not been tried yet.

b. Compute b

=(f

F

−1

)B.

c. Compute H =J ∧unit(b

d. If (H≥ε),setJ ←unit(H).

6. Return σ J.

The implementation of this algorithm can be made very efﬁcient by generating

optimized code for each combination of arguments in Steps 5b and 5c, see Sect. 4.

The algorithm constantly forces blades and vectors to unit size in Steps 5c and 5d.

In an optimized version of the algorithm, one may remove these normalizations and

only normalize the ﬁnal outcome. To apply this optimization, modify the following

steps of the algorithm:

Step 5(c). Compute H =J ∧b

Step 5(d). If (

H

Jb



≥ε),setJ ←H.

Step 6. Return σ unit(J).

However, if this optimization is employed, one should be careful that ﬂoating point

precision does not underﬂow, as J may be scaled by a factor of ε b

 in each loop.

3.2 Computational Example

As an illustration, lets us compute the join(A, B) using the FastJoin algorithm using

A =e

∧e

B =e

∧(0.8e

+0.3e

+0.5e

) =0.3e

∧e

−0.8e

∧e

+0.5e

∧e

Note that in this example, values are rounded to two decimals.

Since the input is nontrivial and the grades of A and B are equal, Steps 1 and 2

of the algorithm have no effect. We continue with Step 3.

Step 3: J ←A.

Step 4: F =−0.8e

∧e

Step 5a, loop 1: f

Step 5b, loop 1: b

=1.00 ∗e

+0.38 ∗e

+0.63 ∗e

Step 5c, loop 1: H ←0.51 ∗e

∧e

+0.30 ∗e

∧e

Step 5d, loop 1: J ←0.86 ∗e

∧e

+0.51 ∗e

∧e

464 D. Fontijne and L. Dorst

Step 5a, loop 2: f

Step 5b, loop 2: b

Step 5c, loop 2: H ←0.

Step 5d, loop 2: H < eps,soJ is not overwritten.

Hence the returned result is the grade 3 blade 0.86∗e

∧e

+0.51∗e

∧e

3.3 Grade Stability of Fast Join Algorithm

The main loop of the above join algorithm tries to increase the grade of the current J

by taking the outer product with factors of B in Step 5c and accepts the result if its

norm is above some threshold value ε. This leads to a problem in stability of the

grade of result: the factorization of B is dependent on the basis, and hence so is the

grade of J as computed by the algorithm.

This is best illustrated using a simple example. Let A and B be the blades

A =e

∧e

, B =e

∧e

+0.2 e

∧e

+0.2 e

∧e

and compute the factors

−0.2 e

, b

+0.2 e

If we now compute Step 5c of the FastJoin algorithm for both b

and b

, we ﬁnd

=A ∧unit(b

) =−0.196e

∧e

=A ∧unit(b

) =0.196e

∧e

Therefore, H

=H

=0.196. If ε is taken as 0.2, the algorithm would accept

neither H

nor H

, and hence return a join of grade 2.

Now suppose that we rotate the input over π/4inthee

∧e

-plane, producing



∧e

, B



∧e

+0.283 e

∧e

This leads to a different factorization,



−0.283 e

, b



resulting in



∧unit







≈−0.272 e

∧e

Hence, with ε =0.2, the algorithm would now return a J of grade 3. Only the posi-

tion of the blades relative to the basis changed relative to the previous situation, not

their mutual position, so we would have expected the grade of their geometrical

join

to have been the same in both cases.

In some application where speed is essential and precision matters less, this be-

havior may be acceptable, but in general one would prefer the grade of the computed

join to be independent of the choice of basis. Hence modiﬁcations to the FastJoin

algorithm are required.

Efﬁcient Algorithms for Factorization and Join of Blades 465

3.4 Improved Fast Join Algorithm

If the required grade of the join J could somehow be computed in advance, this

knowledge could be used to guide the FastJoin algorithm. For instance, we could

ﬁrst run the FastJoin algorithm as is, and if the resulting J does not have the required

grade, we could “lower our standards” and selectively accept some b

that were

rejected earlier, until the grade of the

join is correct. To compute the required grade

of the

join, we use the following equation [6]:

grade



join(A, B)



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

. (4)

The delta product Δ is the geometric symmetric difference [4]. It can be computed

as the highest-grade part of the geometric product AB which is nonzero. In prac-

tice, however, we have to deal with ﬂoating point round-off errors. To quantify this,

when selecting the top-grade part of the geometric product, we use a small thresh-

old δ and select as the value of AΔB the highest-grade part of (unit(A)unit(B))

which has a norm ≥δ. The grade of the delta product can be determined efﬁciently

using lazy evaluation; details follow in Sect. 4.4, for now we call this function

FastDeltaGrade(A, B,δ).

We use this function to improve our FastJoin algorithm as follows:

Algorithm StableFastJoin(A, B,ε,δ):

Start with Steps 1–5 of FastJoin(A, B,ε).

6. If (grade(J) =n) or (grade(J) =grade(A) +grade(B)), return σ J.(σ was com-

puted in Step 2 of the FastJoin algorithm.) Otherwise:

7. Compute grade(

join(A, B)) using (4) and threshold δ for the delta product.

8. While (grade(J)<grade(

join(A, B)))

a. For all valid i, compute b

=(f

F

−1

)B.

Set b

to that b

which leads to the largest J ∧b

.

b. Update J ←J ∧b

9. Return σ J.

Note that one should set ε ≥δ, or else in Step 7 blade J may already have a grade

which is larger than the grade required by FastDeltaGrade(A, B,δ). Also note that

in Step 8a it makes no sense to try any of the b

which were already accepted in

Step 5d or in an earlier iteration of Step 8a.

3.5 Numerical Stability of the Fast Join Algorithms

Another issue in the FastJoin algorithm is numerical stability due to the use of ﬂoat-

ing point values which causes round off errors.

466 D. Fontijne and L. Dorst

The main thing we can do to improve numerical stability of both join algorithms

is to choose good factors b

in Step 5c of the algorithms. Each factor b

potentially

extending the current iteration value of the

join J can be written as the sum



⊥

where b



is parallel to the current J, and b

⊥

is orthogonal to it. The outer product

rejects b



and uses only b

⊥

, i.e.,

H =J ∧b

=J ∧b

⊥

However, when b



b

⊥

, it is likely that the ﬂoating point precision of J ∧b

is low. Hence we would prefer to select those b

which result in the largest J ∧b



because that b

is most orthogonal to the current J.

One solution is trying every b

in each loop and use the one which results in the

largest norm, but this is inefﬁcient.

Fortunately, it is trivial to adjust the StableFastJoin algorithm to be more pre-

cise, with only a minimal performance impact. By running Steps 1 through 5 (its

FastJoin part) using a relatively large threshold like ε = 10

−2

, we only accept fac-

tors in Step 5d which are reasonably orthogonal to J. When the input blades are

in a nondegenerate conﬁguration, the computed J will have the required grade (as

veriﬁed by Step 6 or 7). Otherwise, we will ﬁnd the best factors in Step 8. This step

is more expensive, but it rarely needs to execute.

It is important to realize that in the StableFastJoin algorithm, δ only affects the

correctness of the grade of result, while ε has an effect on the efﬁciency of com-

putation and the numerical precision of the result. That is, a larger ε will let the

algorithm run more slowly on average, but it will also increase numerical precision

of the outcome.

4 Implementation

Below we describe the interesting parts of the implementation of each algorithm.

First, we motivate why code generation was used to write the implementations. Then

we describe the implementation in some more detail, also showing a small piece of

code from the core of each, in order to convey the essence of each operation at the

level of coordinates. Because these functions are used internally in the implementa-

tion, coordinates are passed in the form of arrays of ﬂoating point values, instead of

being encapsulated in classes as is usual in Gaigen 2.

We also give a complexity order of the size of the generated code. This is relevant

because the amount of generated code is large enough to limit the practicality of the

code generation approach, especially for the

join algorithms.

Efﬁcient Algorithms for Factorization and Join of Blades 467

4.1 Code Generation

After the FastFactorization algorithm has found the largest basis blade F, it only

has to scale, copy, and selectively negate certain coordinates from the input blade

to the output factors. Exactly what coordinates have to be copied where depends

on the grade of the input blade and which basis blade was the largest. This can be

implemented very efﬁciently by explicitly spelling out the code for each possible

case. However, writing that amount of code by hand is tedious and error-prone.

Thus, to implement the FastFactorization and the FastJoin algorithms, we have

written a code generator on top of the Gaigen 2 code generation framework [10].

The code generator generates a C++ implementation of the algorithms for a speciﬁc

n and for a speciﬁc order and orientation of the basis elements. It generates the code

of the factorization code without any conditionals (if, else, switch) for each

valid case. For the

join, it also writes out the outer products. We have also written a

code generator for evaluation of the grade of the delta product.

4.2 Implementation of the Fast Factorization Algorithm

The generated implementation of the FastFactorization algorithm ﬁrst ﬁlters out

trivial cases (0, scalars, vectors, pseudoscalars). These are handled separately.

If the input blade is not one of these trivial cases, the implementation ﬁnds the

absolute largest coordinate of the input blade and the respective basis blade F.The

proper scale β of the factorization is computed, and the coordinates are “normal-

ized,” so that the absolute largest coordinate is 1.

The implementation then calls an optimized factorization function which imple-

ments the actual factorization. One such generated function is available for each

possible largest basis blade F, and they are called via a lookup table. Figure 2 shows

an example of such a function. It is clearly visible that the function just copies (and

possibly negates) coordinates of the input blade to the coordinates of a factor. Also,

a number of coordinates of the factors are set to one or zero, as was already apparent

in (3).

The delightful computational simplicity of this way of factorization is the main

idea of this paper. It makes our factorization fast and, with it, the algorithms for the

join that employ it.

In the order of







of these functions are generated, each with code size

proportional to kn, for a total code size of O(







) =O(n

n−1

4.3 Implementation of the Fast Join Algorithm

The generated implementation of the FastJoin algorithm closely follows its descrip-

tion in Sect. 3.1.