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.
