diff options
| -rw-r--r-- | R/optimization.R | 217 |
1 files changed, 107 insertions, 110 deletions
diff --git a/R/optimization.R b/R/optimization.R index 8f662b4e..b05a9cec 100644 --- a/R/optimization.R +++ b/R/optimization.R @@ -1,5 +1,4 @@ #various optimization functions
-library(Rglpk)
"blkdiag" <-
function (m1, m2, p.tr = 0, p.ll = 0)
@@ -14,6 +13,7 @@ function (m1, m2, p.tr = 0, p.ll = 0) }
linfit<-function(A, q, gamma){
+ require(Rglpk)
#minimize ||Ax-q||_1 + gamma||DAx||_1
#1^Tx=1
#x>=0
@@ -40,7 +40,7 @@ linfit1<-function(P, b, A, q, gamma){ # s.t. Px = b
#1^Tx=1
#x>=0
-
+ require(Rglpk)
n<-nrow(P)
p<-ncol(P)
r<-length(q)
@@ -64,7 +64,7 @@ linfit2<-function(P, A, q, gamma){ # s.t. Px=0
#1^Tx=1
#x>=0
-
+ require(Rglpk)
n <- nrow(P)
p <- ncol(P)
b <- rep(0,n)
@@ -305,134 +305,131 @@ KLfit2<-function(P,q,b=rep(0,nrow(P))){ }
solveKKT<-function(H,A,g,h,inv){
- #solves the KKT system
- #inv(H,v) should solve H^-1g
- S <- -t(A)%*%inv(H,A)
- if(ncol(A)==1){
- Dy <- 1/S*(h-t(A)%*%inv(H,g))
- }else{
- Dy <- solve(S,h-t(A)%*%inv(H,g))
- }
- Dx <- inv(H,g-A%*%Dy)
- return(list(Dx=as.vector(Dx),Dy=as.vector(Dy)))
+ ## solves the KKT system
+ ## inv(H,v) should solve H^-1g
+ S <- -t(A)%*%inv(H,A)
+ if(ncol(A)==1){
+ Dy <- 1/S*(h-t(A)%*%inv(H,g))
+ }else{
+ Dy <- solve(S,h-t(A)%*%inv(H,g))
+ }
+ Dx <- inv(H,g-A%*%Dy)
+ return(list(Dx=as.vector(Dx),Dy=as.vector(Dy)))
}
lmconst<-function(Y,X,A,b,pos){
- #solves for \beta min ||Y-X\beta||
- # s.t. A\beta = b
- # \beta>=0 (if pos = TRUE)
- alpha<-0.4
- beta<-0.8
- flag<-1
- #init x and nu
- x<-rep(1,ncol(A))
- nu<-rep(1,nrow(A))
+ ##solves for \beta min ||Y-X\beta||
+ ## s.t. A\beta = b
+ ## \beta>=0 (if pos = TRUE)
+ alpha <- 0.4
+ beta <- 0.8
+ flag <- 1
+ ##init x and nu
+ x <- rep(1,ncol(A))
+ nu <- rep(1,nrow(A))
- t<-pos
- XTX<-crossprod(X,X)
+ t <- pos
+ XTX <- crossprod(X,X)
while(t >1e-10 || flag){
- t <- alpha*t
- flag <- 0
- decr<-Inf
- niter<-1
- while(decr>1e-6 & niter<100){
- #solve the KKT equation
- rdual<--crossprod(X,Y-X%*%x)+crossprod(A,nu)-t/x
- rprimal<-A%*%x-b
- H<-XTX+t*diag(1/x^2)
- KKT<-solveKKT(H,t(A),-rdual,-rprimal,invchol)
- #KKT<-rbind(cbind(H,t(A)),cbind(A,matrix(0,nrow(A),nrow(A))))
- #temp<--solve(KKT,rbind(rdual,rprimal))
- #Dx_pd<-temp[1:nrow(H)]
- #Dnu_pd<-temp[-(1:nrow(H))]
- Dnu_pd<-KKT$Dy
- Dx_pd<-KKT$Dx
- #backtracking search
- phi<-1
- newx<-x+phi*Dx_pd
- newnu<-nu+phi*Dnu_pd
- if( pos==1 ){
- while(min(newx)<0){
- phi<-phi*beta
- newx<-x+phi*Dx_pd
- newnu<-nu+phi*Dnu_pd
- }
- }
- newrdual<--crossprod(X,Y-X%*%newx)+crossprod(A,newnu)-t/newx
- newrprimal<-A%*%newx-b
- while((sqrt(sum(newrdual^2)+sum(newrprimal^2))>(1-alpha*phi)*sqrt(sum(rdual^2)+sum(rprimal^2)))){
- phi<-phi*beta
- newx<-x+phi*Dx_pd
- newnu<-nu+phi*Dnu_pd
- newrdual<--crossprod(X,Y-X%*%newx)+crossprod(A,newnu)-t/newx
- newrprimal<-A%*%newx-b
- }
- x<-newx
- nu<-newnu
- decr<-sqrt(sum(newrdual^2)+sum(newrprimal^2))
- #cat(decr,"\n")
- niter<-niter+1
- }
+ t <- alpha*t
+ flag <- 0
+ decr<-Inf
+ niter<-1
+ while(decr>1e-6 & niter<100){
+ ## solve the KKT equation
+ rdual <- -crossprod(X,Y-X%*%x) + crossprod(A,nu) - t/x
+ rprimal <- A%*%x-b
+ H <- XTX+t*diag(1/x^2)
+ KKT <- solveKKT(H,t(A),-rdual,-rprimal,invchol)
+ Dnu_pd <- KKT$Dy
+ Dx_pd <- KKT$Dx
+ ## backtracking search
+ phi<-1
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ if( pos==1 ){
+ while(min(newx)<0){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ }
+ }
+ newrdual<--crossprod(X,Y-X%*%newx)+crossprod(A,newnu)-t/newx
+ newrprimal<-A%*%newx-b
+ while((sqrt(sum(newrdual^2)+sum(newrprimal^2))>(1-alpha*phi)*sqrt(sum(rdual^2)+sum(rprimal^2)))){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ newrdual<--crossprod(X,Y-X%*%newx)+crossprod(A,newnu)-t/newx
+ newrprimal<-A%*%newx-b
+ }
+ x<-newx
+ nu<-newnu
+ decr<-sqrt(sum(newrdual^2)+sum(newrprimal^2))
+ ##cat(decr,"\n")
+ niter<-niter+1
+ }
}
- return(list(obj=norm(Y-X%*%x),beta=x))
+ return(list(obj=norm(Y-X%*%x),beta=x))
}
#cleanup cdf
cleancdf<-function(x,lambda){
- #minimize for a ||a-x||_1+lambda*||diff(x)||_1
- #subject to diff(x)>=0
- r<-length(x)
- obj<-c(rep(0,r),rep(1,r), rep(lambda,r-1))
- D<-diff(diag(r))
- const.mat<-rbind(cbind(diag(r),diag(r),matrix(0,r,r-1)),
- cbind(D,matrix(0,r-1,r),matrix(0,r-1,r-1)),
- cbind(D,matrix(0,r-1,r),-diag(r-1)),
- cbind(diag(r),-diag(r),matrix(0,r,r-1) ))
+ require(Rglpk)
+ ##minimize for a ||a-x||_1+lambda*||diff(x)||_1
+ ##subject to diff(x)>=0
+ r<-length(x)
+ obj<-c(rep(0,r),rep(1,r), rep(lambda,r-1))
+ D<-diff(diag(r))
+ const.mat<-rbind(cbind(diag(r),diag(r),matrix(0,r,r-1)),
+ cbind(D,matrix(0,r-1,r),matrix(0,r-1,r-1)),
+ cbind(D,matrix(0,r-1,r),-diag(r-1)),
+ cbind(diag(r),-diag(r),matrix(0,r,r-1) ))
- const.dir<-c(rep(">=",2*r-1),rep("<=",2*r-1))#"
- const.rhs<-c(x,rep(0,r-1),rep(0,r-1),x)
- model <- Rglpk_solve_LP(obj,const.mat,const.dir,const.rhs)
- model$solution[1:r]
+ const.dir<-c(rep(">=",2*r-1),rep("<=",2*r-1))#"
+ const.rhs<-c(x,rep(0,r-1),rep(0,r-1),x)
+ model <- Rglpk_solve_LP(obj,const.mat,const.dir,const.rhs)
+ model$solution[1:r]
}
posls<-function(Y,X,w=rep(1,nrow(X))){
- #minimize ||Y-X\beta||_w
- #subject to beta>=0
- beta <- rep(1,ncol(X))
- t <- 1
- phi <- 0.6
- for(i in 1:30){
- dbeta<-rep(Inf,ncol(X))
- while(norm(dbeta)>=1e-10){
- g<--crossprod(X,w*(Y-X%*%beta))-t/beta
- H<-crossprod(X,sweep(X,1,w,"*"))+diag(as.vector(t/beta^2))
- dbeta<-solve(H,g)
- while(min(beta-dbeta)<0){
- dbeta<-phi*dbeta
- }
- beta<-beta-dbeta
- cat(norm(dbeta),"\n")
+ ##minimize ||Y-X\beta||_w
+ ##subject to beta>=0
+ beta <- rep(1,ncol(X))
+ t <- 1
+ phi <- 0.6
+ for(i in 1:30){
+ dbeta<-rep(Inf,ncol(X))
+ while(norm(dbeta)>=1e-10){
+ g<--crossprod(X,w*(Y-X%*%beta))-t/beta
+ H<-crossprod(X,sweep(X,1,w,"*"))+diag(as.vector(t/beta^2))
+ dbeta<-solve(H,g)
+ while(min(beta-dbeta)<0){
+ dbeta<-phi*dbeta
+ }
+ beta<-beta-dbeta
+ cat(norm(dbeta),"\n")
}
- t<-phi*t
- }
- return( beta )
+ t<-phi*t
+ }
+ return( beta )
}
invchol <- function(H,b){
- R <- chol(H)
- z <- forwardsolve(t(R),b)
- x <- backsolve(R,z)
- return(x)
+ R <- chol(H)
+ z <- forwardsolve(t(R),b)
+ x <- backsolve(R,z)
+ return(x)
}
-#solve the KKT equation by block elemination
+#solve the KKT equation by block elimination
solveswm <- function(A,B,C,b){
- #solve Dx=b for D=A+BC with A diagonal and B and C have low ranks
- z <- 1/A*b
- E <- diag(nrow(C))+C%*%sweep(B,1,A,"/")#"
- w <- solve(E,C%*%z)
- x <- z-sweep(B%*%w,1,A,"/")
- return( x )
+ ##solve Dx=b for D=A+BC with A diagonal and B and C have low ranks
+ z <- 1/A*b
+ E <- diag(nrow(C))+C%*%sweep(B,1,A,"/")#"
+ w <- solve(E,C%*%z)
+ x <- z-sweep(B%*%w,1,A,"/")
+ return( x )
}
norm <- function(x){
|