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!
>
> 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()? ***
>
After having a closer look at this problem, I believe you did not include the
constraint x_i >= 0 in the call to solve.QP.
So with this modification of your code
A <- matrix(rep(1,3),nrow=4,ncol=3,byrow=TRUE)
A[2:4,] <- diag(3)
b <- c(1,0,0,0)
sol3 <- solve.QP(D, d, t(A), b, meq = 1) # first row of A is an equality
sol3
p0 <- c(C %*% sol3$solution)
r0 <- sqrt(-sol3$value)
p0
r0
sqrt(colSums((C - p0)^2))
one gets the correct answer.
BTW LowRankQP seems to postulate x_i >=0 if I read its manual correctly.
Berend
______________________________________________
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.
