I wanted to solve the following geometric optimization problem, sometimes
called the "enclosing ball problem":
Given a set P = {p_1, ..., p_n} of n points in R^d, find a point p_0 such
that max ||p_i - p_0|| is minimized.
A known algorithm to solve this as a Qudratic Programming task is as follows:
Define matrix C = (p_1, ..., p_n), i.e. coordinates of points in columns,
and minimize
x' C' C x - \sum (p_i' p_i) x_i
subject to
\sum x_1 = 1, x_i >= 0 .
Let x0 = (x0_1, ..., x0_n) be an optimal solution, then the linear combination
p0 = \sum x0_i p_i
is the center of the smallest enclosing ball, and the negative value of the
objective function at x0 is the squared radius of the ball.
As an example, find the smallest ball enclosing the points (1,0,0), (0,1,0),
and (0.5, 0.5, 0.1) in 3-dim. space. We would expect the center of the ball to
be at (0.5, 0.5, 0), and with radius 1/sqrt(2).
C <- matrix( c(1, 0, 0.5,
0, 1, 0.5,
0, 0, 0.1), ncol = 3, byrow = TRUE)
# We need to define D = 2 C' C, d = (p_i' p_i), A = c(1,1,1), and b = 1
D <- 2 * t(C) %*% C
d <- apply(C^2, 2, sum)
A <- matrix(1, nrow = 1, ncol = 3)
b <- 1
# Now solve with solve.QP() in package quadprog ...
require(quadprog)
sol1 <- solve.QP(D, d, t(A), b, meq = 1)
# ... and check the result
sum(sol1$solution) # 1
p0 <- c(C %*% sol1$solution) # 0.50 0.50 -2.45
r0 <- sqrt(-sol1$value) # 2.55
# distance of all points to the center
sqrt(colSums((C - p0)^2)) # 2.55 2.55 2.55
As a result we get a ball such that all points lie on the boundary.
The same happens for 100 points in 100-dim. space (to keep D positive
definite, n = d is required).
That is a very nice, even astounding result, but not what I had hoped for!
Compare this with another quadratic programming solver in R, for example
LowRankQP() in the package of the same name.
require(LowRankQP)
sol2 <- LowRankQP(D, -d, A, b, u = rep(3, 3), method="LU")
p2 <- c(C %*% sol2$alpha) # 5.000000e-01 5.000000e-01 1.032019e-12
sqrt(colSums((C - p2)^2)) # 0.7071068 0.7071068 0.1000000
The correct result, and we get correct solutions in higher dimensions, too.
LowRankQP works also if we apply it with n > d, e.g., 100 points in 2-
or 3-dimensional space, when matrix D is not positive definite -- solve.QP( )
does not work in these cases.
*** What do I do wrong in calling solve.QP()? ***
My apologies for this overly long request. I am sending it through the weekend
hoping that someone may have a bit of time to look at it more closely.
Hans Werner
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
