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 > ______________________________________________ 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.