Fast Smallest-Enclosing-Ball Computation in High Dimensions Martin - - PowerPoint PPT Presentation

▶
fast smallest enclosing ball computation in high
SMART_READER_LITE
LIVE PREVIEW

Fast Smallest-Enclosing-Ball Computation in High Dimensions Martin - - PowerPoint PPT Presentation

Feb 14, 2024 89 likes •827 views

Fast Smallest-Enclosing-Ball Computation in High Dimensions Martin Kutz FU Berlin joint work with Bernd G artner and Kaspar Fischer, ETH Z urich Martin Kutz Smallest Enclosing Balls The Smallest Enclosing Ball given: finite point set

slide-1
SLIDE 1

Fast Smallest-Enclosing-Ball Computation in High Dimensions

Martin Kutz

FU Berlin joint work with Bernd G¨ artner and Kaspar Fischer, ETH Z¨ urich

Martin Kutz — Smallest Enclosing Balls

slide-2
SLIDE 2

The Smallest Enclosing Ball

given: finite point set S in Rd

Martin Kutz — Smallest Enclosing Balls 1

slide-3
SLIDE 3

The Smallest Enclosing Ball

given: finite point set S in Rd wanted: smallest ball B = B(c, r) := {x : ||x − c|| ≤ r} containing S

Martin Kutz — Smallest Enclosing Balls 2

slide-4
SLIDE 4

The Smallest Enclosing Ball

given: finite point set S in Rd wanted: smallest ball B = B(c, r) := {x : ||x − c|| ≤ r} containing S Call this unique minimal B the smallest enclosing ball of S, denoted seb(S).

Martin Kutz — Smallest Enclosing Balls 3

slide-5
SLIDE 5

Applications

  • visibility culling and bounding sphere hierarchies in 3D computer graphics
  • clustering (e.g. for support-vector machines) — many dimensions
  • nearest neighbor search

Martin Kutz — Smallest Enclosing Balls 4

slide-6
SLIDE 6

Previous Work

  • Welzl proposed randomized combinatorial algorithm, implemented by

G¨ artner, fast for d ≤ 30, impractical above.

  • Quadratic-programming approach by G¨

artner & Sch¨

  • nherr, uses exact

arithmetic, up to d = 300.

Martin Kutz — Smallest Enclosing Balls 5

slide-7
SLIDE 7

Previous Work

  • Welzl proposed randomized combinatorial algorithm, implemented by

G¨ artner, fast for d ≤ 30, impractical above.

  • Quadratic-programming approach by G¨

artner & Sch¨

  • nherr, uses exact

arithmetic, up to d = 300.

  • General-purpose QP-solver CPLEX, solves d ≤ 3,000.

Martin Kutz — Smallest Enclosing Balls 6

slide-8
SLIDE 8

Previous Work

  • Welzl proposed randomized combinatorial algorithm, implemented by

G¨ artner, fast for d ≤ 30, impractical above.

  • Quadratic-programming approach by G¨

artner & Sch¨

  • nherr, uses exact

arithmetic, up to d = 300.

  • General-purpose QP-solver CPLEX, solves d ≤ 3,000.
  • Zhou, Toh, and Sun use interior-point method to find approximate

solution, up to d = 10,000.

  • Kumar, Mitchell, Yildrum compute approximation with core sets, results

given up to d = 1,400.

Martin Kutz — Smallest Enclosing Balls 7

slide-9
SLIDE 9

Our Algorithm

  • simple combinatorial algorithm

(not approximation)

Martin Kutz — Smallest Enclosing Balls 8

slide-10
SLIDE 10

Our Algorithm

  • simple combinatorial algorithm

(not approximation)

  • similar to LP simplex-method
  • equipped with pivot scheme to avoid cycling

Martin Kutz — Smallest Enclosing Balls 9

slide-11
SLIDE 11

Our Algorithm

  • simple combinatorial algorithm

(not approximation)

  • similar to LP simplex-method
  • equipped with pivot scheme to avoid cycling
  • C++ floating-point implementation: solves several thousand points in a

few thousand dimensions

Martin Kutz — Smallest Enclosing Balls 10

slide-12
SLIDE 12

Our Algorithm

  • simple combinatorial algorithm

(not approximation)

  • similar to LP simplex-method
  • equipped with pivot scheme to avoid cycling
  • C++ floating-point implementation: solves several thousand points in a

few thousand dimensions

  • idea not completely new; Hopp & Reeve presented similar algorithm but

without proofs, some details unclear, 3D only

Martin Kutz — Smallest Enclosing Balls 11

slide-13
SLIDE 13

The Basic Idea: Deflating a Ball

Iteratively shrink an enclosing Ball B = B(c, T) represented by

  • a current center c,
  • an affinely independent subset T ⊆ S of points at a common distance

from c — the support set

Martin Kutz — Smallest Enclosing Balls 12

slide-14
SLIDE 14

The Basic Idea: Deflating a Ball

Iteratively shrink an enclosing Ball B = B(c, T) represented by

  • a current center c,
  • an affinely independent subset T ⊆ S of points at a common distance

from c — the support set Invariants: T ⊂ ∂B(c, T) S ⊂ B(c, T)

Martin Kutz — Smallest Enclosing Balls 13

slide-15
SLIDE 15

2D “Example”

Martin Kutz — Smallest Enclosing Balls 14

slide-16
SLIDE 16

2D “Example”

Martin Kutz — Smallest Enclosing Balls 15

slide-17
SLIDE 17

2D “Example”

Martin Kutz — Smallest Enclosing Balls 16

slide-18
SLIDE 18

2D “Example”

Martin Kutz — Smallest Enclosing Balls 17

slide-19
SLIDE 19

2D “Example”

Martin Kutz — Smallest Enclosing Balls 18

slide-20
SLIDE 20

2D “Example”

Martin Kutz — Smallest Enclosing Balls 19

slide-21
SLIDE 21

2D “Example”

Martin Kutz — Smallest Enclosing Balls 20

slide-22
SLIDE 22

2D “Example”

Martin Kutz — Smallest Enclosing Balls 21

slide-23
SLIDE 23

Termination Criterion

Martin Kutz — Smallest Enclosing Balls 22

slide-24
SLIDE 24

Termination Criterion

Lemma (Seidel). Let T be set of points on boundary

  • f some ball B with center c.

Then B = seb(T) ⇐ ⇒ c ∈ conv(T).

Martin Kutz — Smallest Enclosing Balls 23

slide-25
SLIDE 25

How to Shrink

Moving Step [ Precondition c ∈ aff(T) ] Move c orthogonally towards aff(T), i.e., heading for closest point in aff(T).

aff(T)

Martin Kutz — Smallest Enclosing Balls 24

slide-26
SLIDE 26

How to Shrink

Moving Step [ Precondition c ∈ aff(T) ] Move c orthogonally towards aff(T), i.e., heading for closest point in aff(T). For any fixed point c′ on this path, T stays on sphere around c′. Especially, our target point is the center of the unique sphere through T in aff(T), called circumcenter of T.

aff(T)

Martin Kutz — Smallest Enclosing Balls 25

slide-27
SLIDE 27

How to Shrink

Moving Step [ Precondition c ∈ aff(T) ] Move c orthogonally towards aff(T), i.e., heading for closest point in aff(T). For any fixed point c′ on this path, T stays on sphere around c′. Especially, our target point is the center of the unique sphere through T in aff(T), called circumcenter of T. Stop movement when shrinking boundary hits new point of S, insert it into T;

  • therwise just stop with c in aff(T).

aff(T)

Martin Kutz — Smallest Enclosing Balls 26

slide-28
SLIDE 28

How to Shrink

Dropping Step Necessary if c ∈ aff(T) \ conv(T). Must remove a point from T.

+ − +

Martin Kutz — Smallest Enclosing Balls 27

slide-29
SLIDE 29

How to Shrink

Dropping Step Necessary if c ∈ aff(T) \ conv(T). Must remove a point from T. Pick one with negative coefficient in affine representation c =

  • p∈T

λpp,

  • p∈T

λp = 1.

+ − +

Martin Kutz — Smallest Enclosing Balls 28

slide-30
SLIDE 30

How to Shrink

Dropping Step Necessary if c ∈ aff(T) \ conv(T). Must remove a point from T. Pick one with negative coefficient in affine representation c =

  • p∈T

λpp,

  • p∈T

λp = 1. Afterwards, c lies outside the new aff(T), so it we can move again. The next move will not recollect the dropped point.

+ − +

Martin Kutz — Smallest Enclosing Balls 29

slide-31
SLIDE 31

The Whole Algorithm

c := any point of S; T := {p}, with some p ∈ S at maximal distance from c; while c ∈ conv(T) do [ Invariant: B(c, T) ⊃ S, ∂B(c, T) ⊃ T, and T affinely independent ] if c ∈ aff(T) then drop T-point with negative coefficient in aff. rep. of c; [ Invariant: c ∈ aff(T) ]

Martin Kutz — Smallest Enclosing Balls 30

slide-32
SLIDE 32

The Whole Algorithm

c := any point of S; T := {p}, with some p ∈ S at maximal distance from c; while c ∈ conv(T) do [ Invariant: B(c, T) ⊃ S, ∂B(c, T) ⊃ T, and T affinely independent ] if c ∈ aff(T) then drop T-point with negative coefficient in aff. rep. of c; [ Invariant: c ∈ aff(T) ] move c towards aff(T), stop when boundary hits new point q ∈ S or c reaches aff(T); if point stopped us then T := T ∪ {q}; end while;

Martin Kutz — Smallest Enclosing Balls 31

slide-33
SLIDE 33

Our Example Again

Martin Kutz — Smallest Enclosing Balls 32

slide-34
SLIDE 34

Our Example Again

Martin Kutz — Smallest Enclosing Balls 33

slide-35
SLIDE 35

Our Example Again

Martin Kutz — Smallest Enclosing Balls 34

slide-36
SLIDE 36

Our Example Again

Martin Kutz — Smallest Enclosing Balls 35

slide-37
SLIDE 37

Our Example Again

Martin Kutz — Smallest Enclosing Balls 36

slide-38
SLIDE 38

Our Example Again

+ − +

Martin Kutz — Smallest Enclosing Balls 37

slide-39
SLIDE 39

Our Example Again

+ +

Martin Kutz — Smallest Enclosing Balls 38

slide-40
SLIDE 40

Our Example Again

Martin Kutz — Smallest Enclosing Balls 39

slide-41
SLIDE 41

Our Example Again

+ + +

Martin Kutz — Smallest Enclosing Balls 40

slide-42
SLIDE 42

Martin Kutz — Smallest Enclosing Balls 41

slide-43
SLIDE 43

Martin Kutz — Smallest Enclosing Balls 42

slide-44
SLIDE 44

Martin Kutz — Smallest Enclosing Balls 43

slide-45
SLIDE 45

+ + −

Martin Kutz — Smallest Enclosing Balls 44

slide-46
SLIDE 46

+ +

Martin Kutz — Smallest Enclosing Balls 45

slide-47
SLIDE 47

Correctness & Termination

Correctness “clear” from invariants.

Martin Kutz — Smallest Enclosing Balls 46

slide-48
SLIDE 48

Correctness & Termination

Correctness “clear” from invariants. while c ∈ conv(T) do [ Invariant: B(c, T) ⊃ S, ∂B(c, T) ⊃ T, and T affinely independent ] if c ∈ aff(T) then drop T-point with negative coefficient in aff. rep. of c; [ Invariant: c ∈ aff(T) ] move c towards aff(T), stop when boundary hits new point q ∈ S or c reaches aff(T); if point stopped us then T := T ∪ {q}; end while;

Martin Kutz — Smallest Enclosing Balls 47

slide-49
SLIDE 49

Correctness & Termination

Correctness “clear” from invariants. Termination more complicated.

Martin Kutz — Smallest Enclosing Balls 48

slide-50
SLIDE 50

Correctness & Termination

Correctness “clear” from invariants. Termination more complicated. Proposition. In the non-degenerate case (no affinely dependent subset T ⊆ S lies on a sphere) the algorithm terminates.

Martin Kutz — Smallest Enclosing Balls 49

slide-51
SLIDE 51

Correctness & Termination

Correctness “clear” from invariants. Termination more complicated. Proposition. In the non-degenerate case (no affinely dependent subset T ⊆ S lies on a sphere) the algorithm terminates. Proof:

  • Negative-coefficient rule prevents immediate re-insertion after drop.
  • Radius decreases after dropping step.
  • At least 1 out of d consecutive iterations performs a drop.
  • Set of all possible balls B(c, T) preceding drops is finite.

Martin Kutz — Smallest Enclosing Balls 50

slide-52
SLIDE 52

How to Prevent Cycling

In degenerate cases cycling may occur, i.e., the center c doesn’t move but

  • nly support set T changes — forever.

Martin Kutz — Smallest Enclosing Balls 51

slide-53
SLIDE 53

How to Prevent Cycling

In degenerate cases cycling may occur, i.e., the center c doesn’t move but

  • nly support set T changes — forever.

Solution: pivot rule, similar to Bland’s rule for simplex algorithm. Index the point set S in arbitrary order. When dropping a point with negative coefficient, pick the one with smallest index. When movement stopped by several points, also pick the one with smallest index.

Martin Kutz — Smallest Enclosing Balls 52

slide-54
SLIDE 54

How to Prevent Cycling

In degenerate cases cycling may occur, i.e., the center c doesn’t move but

  • nly support set T changes — forever.

Solution: pivot rule, similar to Bland’s rule for simplex algorithm. Index the point set S in arbitrary order. When dropping a point with negative coefficient, pick the one with smallest index. When movement stopped by several points, also pick the one with smallest index. Theorem. Using “Bland’s rule” our algorithm terminates.

Martin Kutz — Smallest Enclosing Balls 53

slide-55
SLIDE 55

Technical Details

Data structure for support set T needed that allows requests

  • compute orthogonal projection onto aff(T)

(for walking),

  • compute affine coefficients of point p ∈ aff(T)

(for dropping)

Martin Kutz — Smallest Enclosing Balls 54

slide-56
SLIDE 56

Technical Details

Data structure for support set T needed that allows requests

  • compute orthogonal projection onto aff(T)

(for walking),

  • compute affine coefficients of point p ∈ aff(T)

(for dropping) and updates

  • insert point into T and
  • delete point from T.

Martin Kutz — Smallest Enclosing Balls 55

slide-57
SLIDE 57

Technical Details

a1 a2 b

Let A :=

  • a1

a2 · · · ar

  • (homogenized).

Martin Kutz — Smallest Enclosing Balls 56

slide-58
SLIDE 58

Technical Details

a1 a2 Ax∗ b

Let A :=

  • a1

a2 · · · ar

  • (homogenized).

Let x∗ ∈ a1, a2, . . . , ak minimize the risidual ||Ax − b||2, then Ax∗ is the orthogonal projection

  • f b onto a1, a2, . . . , ak.

Martin Kutz — Smallest Enclosing Balls 57

slide-59
SLIDE 59

Technical Details

a1 a2 Ax∗ b

Let A :=

  • a1

a2 · · · ar

  • (homogenized).

Let x∗ ∈ a1, a2, . . . , ak minimize the risidual ||Ax − b||2, then Ax∗ is the orthogonal projection

  • f b onto a1, a2, . . . , ak.

If b ∈ a1, a2, . . . , ak then the coefficients of x∗ simply are the coefficients of b.

Martin Kutz — Smallest Enclosing Balls 58

slide-60
SLIDE 60

Technical Details

We compute the x∗ that minimizes ||Ax − b|| with QR-decomposition

= R × Q A

Martin Kutz — Smallest Enclosing Balls 59

slide-61
SLIDE 61

Technical Details

We compute the x∗ that minimizes ||Ax − b|| with QR-decomposition

= R × Q A

“Solve” Ax = b via QRx = b ⇐ ⇒ Rx = QTb.

Martin Kutz — Smallest Enclosing Balls 60

slide-62
SLIDE 62

Technical Details

We compute the x∗ that minimizes ||Ax − b|| with QR-decomposition

= R × Q A

“Solve” Ax = b via QRx = b ⇐ ⇒ Rx = QTb. Let y :≈ QTb with last entries zeroed and then solve Rx = y via back substitution.

Martin Kutz — Smallest Enclosing Balls 61

slide-63
SLIDE 63

Technical Details

Update QR-decomposition via Givens rotations.

= R × Q A

Martin Kutz — Smallest Enclosing Balls 62

slide-64
SLIDE 64

Technical Details

Update QR-decomposition via Givens rotations.

= R × Q A

  • add new column a to A (and rotated column QTa to R)

Martin Kutz — Smallest Enclosing Balls 63

slide-65
SLIDE 65

Technical Details

Update QR-decomposition via Givens rotations.

= R × Q A

  • add new column a to A (and rotated column QTa to R)
  • rotate rows of R and colums of Q to reduce R to triangular shape

Martin Kutz — Smallest Enclosing Balls 64

slide-66
SLIDE 66

Technical Details

Update QR-decomposition via Givens rotations.

= R × Q A

  • add new column a to A (and rotated column QTa to R)
  • rotate rows of R and colums of Q to reduce R to triangular shape

Martin Kutz — Smallest Enclosing Balls 65

slide-67
SLIDE 67

Our Implementation

  • single iteration in O(nd) time

Martin Kutz — Smallest Enclosing Balls 66

slide-68
SLIDE 68

Our Implementation

  • single iteration in O(nd) time
  • C++ floating-point
  • Bland’s rule replaced by numerically more stable heuristic

Martin Kutz — Smallest Enclosing Balls 67

slide-69
SLIDE 69

Our Implementation

  • single iteration in O(nd) time
  • C++ floating-point
  • Bland’s rule replaced by numerically more stable heuristic
  • QR decomposition numerically very stable
  • very accurate results, about 1,000 times machine precision

Martin Kutz — Smallest Enclosing Balls 68

slide-70
SLIDE 70

Martin Kutz — Smallest Enclosing Balls 69

slide-71
SLIDE 71

Kumar et al.: n = 1000 Zhou et al.: n = 1000

  • ur algorithm: n = 1000
  • ur algorithm: n = 2000

Uniform Distribution dimension seconds 2000 1500 1000 500 400 350 300 250 200 150 100 50

Martin Kutz — Smallest Enclosing Balls 70

slide-72
SLIDE 72

CPLEX: n = 1000 CPLEX: n = 2000

  • ur algorithm: n = 1000
  • ur algorithm: n = 2000

Almost Spherical Distribution dimension seconds 2000 1500 1000 500 3500 3000 2500 2000 1500 1000 500

Martin Kutz — Smallest Enclosing Balls 71