[R] Optimization of constrained linear least-squares problem

[Stefaan Lhermitte]

Thanx Dimitris, Patrick and Berwin!

For other people interested in this problem, here are two valid 
solutions that work.

1) Re-parameterize  e.g.,

    EM <- c(100,0,0,0,100,0,0,0,100)
    W <- array(EM, c(3,3))
    d <- c(10, 20, 70)

    fn <- function(x){
          x <- exp(x) / sum(exp(x))
          r <- W%*%x - d
        crossprod(r, r)[1,1]
    }
    opt <- optim(rnorm(3), fn)
    res <- exp(opt$par) / sum(exp(opt$par))
    res

    "The first line of the `fn()' function is just a re-pameterization
    of your problem, i.e., if `y' is a vector of real numbers, then it
    is straightforward to see that `x = exp(y) / sum(exp(y))' will be
    real numbers in (0, 1) for which `sum(y)=1'. So instead of finding
    xs that minimize your function under the constraint (which is more
    difficult) you just find the ys using the above transformation." (I
    owe you a drink Dimitris !!!)

2)  Or minimize it as a quadratic function under a linear constraint:

    EM <- c(100,0,0,0,100,0,0,0,100)
    W <- array(EM, c(3,3))
    d <- c(10, 20, 70)
    library(quadprog)
    Dmat <- crossprod(W,W)
    dvec <- crossprod(d,W)
    A <-matrix(c(1,1,1),ncol=1)
    bvec <- 1
    solve.QP(Dmat, dvec, A, bvec, meq=1)

    This is based on the objective function (i.e. the thing you want to
    minimise) : min x'C'Cx - 2 d'Cx + d'd where sum(x) = 1
    (Thanx Berwin!!)

Kind regards,
Stef