[R] The smallest enclosing ball problem

[Hans W.Borchers]

Berend Hasselman <bhh <at> xs4all.nl> writes:
> Forgot to forward my answer to R-help.
> 
> Berend

Thanks, Berend, for thinking about it. \sum xi = 1  is a necessary condition 
to generate a valid geometric solution. The three points in the example are 
very regular and your apporach works, but imagine some random points:

    set.seed(8237)
    C <- matrix(runif(9), 3, 3)
    D <- 2 * t(C) %*% C
    d <- apply(C^2, 2, sum)
    A <- diag(3)
    b <- rep(0,3)

    require(quadprog)
    sol1 <- solve.QP(D, d, A, b, meq = 0)
    sol1                                # now \sum xi = 1is not fulfilled

    p0 <- c(C %*% sol1$solution)        # 0.3707410 0.3305265 0.2352084
    r0 <- sqrt(-sol1$value)             # 0.5495631

    # distance of all points to the center
    sqrt(colSums((C - p0)^2))           # 0.5495631 0.3119314 0.5495631

Unfortunately, this is not the smallest enclosing ball.
LowRankQP will find the true solution with radius 0.3736386 !

    require(LowRankQP)
    A <- matrix(1, nrow = 1, ncol = 3)
    b <- 1
    
    sol2 <- LowRankQP(D, -d, A, b, u = rep(1, 3), method="LU")
  
    p2 <- c(C %*% sol2$alpha)   # 0.5783628 0.5372570 0.2017087
    sqrt(colSums((C - p2)^2))   # 0.3736386 0.3736386 0.3736386

But the strangest thing is that with \sum xi =1 solve.QP positions all points
on the boundary, though (in my opinion) no constraint requests it. So the
question remains:
                  *** What do I do wrong in calling solve.QP()? ***

Hans Werner

> At the risk of making a complete fool of myself.
> 
> Why the restriction:  \sum x_1 = 1, x_i >= 0 ?
> 
> Isn’t just  x_i >= 0 sufficient?
> 
> Change your code with this
> 
> A <- diag(3)
> b <- rep(0,3)
> sol2 <- solve.QP(D, d, A, b, meq = 0)
> sol2
> 
> resulting in this
> 
> $solution
> [1] 0.5 0.5 0.0
> $value
> [1] -0.5
> $unconstrained.solution
> [1]  12.75  12.75 -24.50
> $iterations
> [1] 2 0
> $Lagrangian
> [1] 0.00 0.00 0.49
> $iact
> [1] 3
> 
> And $solution seems to be what you want.
> 
> And:
> 
> p0 <- c(C %*% sol2$solution)
> r0 <- sqrt(-sol2$value)
> 
> # distance of all points to the center
> sqrt(colSums((C - p0)^2))
> 
> gives the same results as LowRankQP for the last expression.
> 
> Berend
>