On 17-11-2013, at 11:32, Hans W.Borchers <hwborch...@googlemail.com> wrote:
> 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 See my second reply to your original post. Modify your code with A <- matrix(rep(1,3),nrow=4,ncol=3,byrow=TRUE) A[2:4,] <- diag(3) b <- c(1,0,0,0) to include constraints x_i>=0 (which LowRankQP includes automatically!) and run solve.QP as follows sol2 <- solve.QP(D, d, t(A), b, meq = 1) sol2 p0 <- c(C %*% sol2$solution) r0 <- sqrt(-sol2$value) p0 r0 # distance of all points to the center sqrt(colSums((C - p0)^2)) and the answers now agree with LowRankQP. 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.
