[R] The smallest enclosing ball problem

[Berend Hasselman]

On 17-11-2013, at 11:32, Hans W.Borchers <hwborchers at googlemail.com> wrote:

> 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

See my second reply to your original post.
Modify your code with

A <- matrix(rep(1,3),nrow=4,ncol=3,byrow=TRUE)
A[2:4,] <- diag(3)   
b <- c(1,0,0,0)

to include constraints x_i>=0 (which LowRankQP includes automatically!) and run solve.QP as follows

sol2 <- solve.QP(D, d, t(A), b, meq = 1)
sol2

p0 <- c(C %*% sol2$solution) 
r0 <- sqrt(-sol2$value) 
p0
r0
# distance of all points to the center
sqrt(colSums((C - p0)^2))
 
and the answers now agree with LowRankQP.

Berend