Forgot to forward my answer to R-help. Berend
On 16-11-2013, at 13:11, Hans W.Borchers <hwborch...@googlemail.com> wrote: > 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! > 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 > 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. ______________________________________________ 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.