The smallest enclosing ball problem

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