Skip to content

Optimization of constrained linear least-squares problem

3 messages · Stefaan Lhermitte, Dimitris Rizopoulos

Original
Stefaan Lhermitte
# 
Dear R-ians,

I want to perform an linear unmixing of image pixels in fractions of 
pure endmembers. Therefore I need to perform a constrained linear 
least-squares problem that looks like :

min || Cx - d || ? where sum(x) = 1.

I have a 3x3 matrix C, containing the values for endmembers and I have a 
3x1 column vector d (for every pixel in the image). In theory my x 
values should all be in the (0,1) interval but I don't want to force 
them so I can check the validity of my solution. I just want to 
calculate the x values. Can anyone help me with this problem? I've been 
checking the optim, optimize, constrOptim and nlm help files, burt  I 
don't understand it very well. Wich function should I use for my 
problem? I did a first test using optim:

# Make my C matrix
EM<- c(4.5000,6.0000,10.5000,5.0000,27.0000,20.7500,16.7500,23.6666,38.7500)
C <- array(EM, c(3,3))

# Take an arbitrary d
d<-c(10, 20, 20)

# Define the function
 fr <- function(x) {
C[1,]*x=d
C[2,]*x=d
C[3,]*x=d
sum(x)=1}

# Perform the optimization
optim(c(0.25,0.25,0.25),fr)

But it did not work. I got the eror couldn't find function. Can anyone 
tell me what functyion I should use for my problem and how should I 
program it?

Thanx in advance,
Stef
Dimitris Rizopoulos
# 
you could re-parameterize, e.g.,

EM <- 
c(4.5000,6.0000,10.5000,5.0000,27.0000,20.7500,16.7500,23.6666,38.7500)
W <- array(EM, c(3,3))
d <- c(10, 20, 20)
##############33
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


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat/
     http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Stefaan Lhermitte" <stefaan.lhermitte at biw.kuleuven.be>
To: <r-help at stat.math.ethz.ch>
Sent: Thursday, March 17, 2005 3:13 PM
Subject: [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