# [R] The smallest enclosing ball problem

Sun Nov 17 11:32:24 CET 2013

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)

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
>  0.5 0.5 0.0
> $value >  -0.5 >$unconstrained.solution
>   12.75  12.75 -24.50
> $iterations >  2 0 >$Lagrangian
>  0.00 0.00 0.49
> $iact >  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
>