# Chapter 2: Introduction to R A <- matrix(c(16,3,2,13,5,10,11,8,9,6,7,12,4,15,14,1), 4,4,byrow=TRUE) A sum(A) colSums(A) rowSums(A) sum(diag(A)) A[1,4]+A[2,3]+A[3,2]+A[4,1] sum(A[outer(1:4,1:4,"+")==5]) B <- rbind(matrix(ceiling(10*runif(10)),2), matrix(rnorm(10),2)) B A <- matrix(c(16,3,2,13,5,10,11,8,9,6,7,12,4,15,14,1), 4,4,byrow=TRUE) A A*A A%*%A t(A)%*%A det(A) eigen(A) U <- eigen(A)$vector L <- diag(eigen(A)$val) U%*%t(U) t(U)%*%U L U%*%L%*%t(U) - A falling <- function(t,gravity=32) gravity*t^2/2 falling(0:5) falling(0:5,8) g <- function(x,a) 0.5*(x + a/x) x <- 1 for(i in 1:10^2) x <- g(x,4) x my.sqrt <- function(a,tol=0/10^6){ x0 <- x <- 1 d <- 10^6 while(d > tol){ x <- g(x0,a) d <- abs(x-x0) x0 <- x } return(x) } my.sqrt(2) sqrt(2) t <- seq(0,2*pi,length=100) y <- sin(t) plot(t,y,type="l") y2 <- sin(t-0.25) lines(t,y2,col=2) y3 <- sin(t-0.5) lines(t,y3,col=3) y <- x <- seq(0,1,by=0.01) z <- outer(x,y,function(x,y) exp(- (x^2 + y^2))) contour(x,y,z) # Chapter 3: Line Search humps <- function(x) 1/((x-0.3)^2+0.01)+1/((x-0.9)^2+0.04)-6 x <- seq(-5,5,by=0.01) plot(x,humps(x),type="l") optimize(humps,c(0.3,1)) abline(v=c(0.3,1),lty=2) humps.plot <- function(x){ y <- 1/((x-0.3)^2 + 0.01) + 1/((x-0.9)^2+0.04) - 6 lines(x,y,col=2,type="h") print(c(x,y)) return(y) } x <- seq(0.3,1,by=0.01) plot(x,humps(x),type="l") abline(v=c(0.3,1),lty=2) optimize(humps.plot,c(0.3,1)) # Chapter 4: Conditions for Unconstraint Solutions Q <- matrix(c(1,-0.5,-0.5,1),2,2) eigen(Q) y <- seq(0,0.5,length=100) x <- seq(0.25,0.75,length=100) z <- outer(x,y,function(x,y) x^2 -x + y + x*y) contour(x,y,z) f.dd <- matrix(c(18,-12,-12,4),2,2) eigen(f.dd) y <- seq(8.5,9.5,length=100) x <- seq(5.5,6.5,length=100) z <- outer(x,y,function(x,y) x^3 -x^2*y + 2*y^2) filled.contour(x,y,z) # Chapter 5: Steepest Decent Q <- matrix(c(0.78,-0.02,-0.12,-0.14,-0.02,0.86,-0.04,0.06, -0.12,-0.04,0.72,-0.08,-0.14,0.06,-0.08,0.74),4,4,byrow=TRUE) Q b <- c(0.76,0.08,1.12,0.68) b quad <- function(x,Q,b) t(x)%*%Q%*%x/2 - sum(x*b) optim(rep(1,4),quad,Q=Q,b=b) solve(Q)%*%b quad.grad <- function(x,Q,b) Q%*%x - b steepest.decent <- function(fun,grad,x,...) { g <- function(a,x,...) fun(x - a*grad(x,...),...) opt <- optimize(g,c(0,10),x=x,...) alpha <- opt$minimum objective <- opt$objective x <- x - alpha*grad(x,...) return(list(x=x,alpha=alpha,objective=objective)) } steepest.decent(quad,quad.grad,rep(1,4),Q=Q,b=b) xn <- x <- rep(1,4) fn <- quad(x,Q,b) for (i in 1:10) { out <- steepest.decent(quad,quad.grad,x,Q=Q,b=b) x <- out$x xn <- cbind(xn,x) fn <- c(fn,out$obj) } fn xn eigen(Q) ((0.94 - 0.52)/(0.94 + 0.52))^2 # Chapter 6: Newton Q <- matrix(c(0.78,-0.02,-0.12,-0.14,-0.02,0.86,-0.04,0.06, -0.12,-0.04,0.72,-0.08,-0.14,0.06,-0.08,0.74), 4,4,byrow=TRUE) b <- c(0.76,0.08,1.12,0.68) quad <- function(x,Q,b) t(x)%*%Q%*%x/2 - sum(x*b) quad.grad <- function(x,Q,b) Q%*%x - b solve(Q)%*%b DFP <- function(fun,grad,x,S,...) { g <- function(a,x,S,...) fun(x - a*S%*%grad(x,...),...) opt <- optimize(g,c(0,10),x=x,S=S,...) alpha <- opt$minimum objective <- opt$objective x1 <- x - alpha*S%*%grad(x,...) dx <- x1-x ddf <- grad(x1,...) - grad(x,...) term1 <- dx %*% t(dx)/sum(dx*ddf) term2 <- S %*% ddf %*% t(ddf) %*% S term2 <- term2/as.numeric(t(ddf) %*% S %*% ddf) S <- S + term1 - term2 x <- x1 return(list(x=x,S=S,alpha=alpha,objective=objective)) } x0 <- rep(1,4) S0 <- diag(4) DFP(quad,quad.grad,x0,S0,Q=Q,b=b) x <- xn <- x0 S <- S0 fn <- quad(x0,Q,b) for (i in 1:3) { out <- DFP(quad,quad.grad,x,S,Q=Q,b=b) x <- out$x S <- out$S xn <- cbind(xn,x) fn <- c(fn,out$obj) } fn xn S solve(Q) # Chapter 8: Quadratic Programming Q.quad <- matrix(c(0.78,-0.02,-0.12,-0.14,-0.02,0.86,-0.04,0.06, -0.12,-0.04,0.72,-0.08,-0.14,0.06,-0.08,0.74), 4,4,byrow=TRUE) Q.quad b.quad <- c(0.76,0.08,1.12,0.68) A.cons <- matrix(c(1,1,1,1,1,1,-1,-1),nrow=2,byrow=TRUE) A.cons b.cons <- c(0,0) b.cons V <- rbind(t(A.cons),matrix(0,nrow(A.cons),nrow(A.cons))) V <- cbind(rbind(Q.quad,A.cons),V) V u <- c(-b.quad,b.cons) u solve(V,u) QP <- function(Q.quad,b.quad,A.cons,b.cons) { d <- length(b.quad) m <- length(b.cons) if(m==1) A.cons <- matrix(A.cons,nrow=1) V <- rbind(t(A.cons),matrix(0,m,m)) V <- cbind(rbind(Q.quad,A.cons),V) u <- c(-b.quad,b.cons) xl <- solve(V,u) x <- xl[1:d] if(m >= 1) lam <- xl[1:m +d] else lam <- NA return(list(x=x,lam=lam)) } QP(Q.quad,b.quad,A.cons,b.cons) V <- rbind(t(A.cons),matrix(0,nrow(A.cons),nrow(A.cons))) V <- cbind(rbind(Q.quad,A.cons),V) V u <- c(-b.quad,b.cons) u solve(V,u) QP <- function(Q.quad,b.quad,A.cons,b.cons) { d <- length(b.quad) m <- length(b.cons) if(m==1) A.cons <- matrix(A.cons,nrow=1) V <- rbind(t(A.cons),matrix(0,m,m)) V <- cbind(rbind(Q.quad,A.cons),V) u <- c(-b.quad,b.cons) xl <- solve(V,u) x <- xl[1:d] if(m >= 1) lam <- xl[1:m +d] else lam <- NA return(list(x=x,lam=lam)) } QP(Q.quad,b.quad,A.cons,b.cons) library(quadprog) ?solve.QP solve.QP(Q.quad,-b.quad,t(A.cons),b.cons,meq=2) solve.QP(Q.quad,-b.quad,t(A.cons),b.cons,meq=2)$solution QP(Q.quad,b.quad,A.cons,b.cons)$x Q.quad <- matrix(c(4,1,1,2),2,2) Q.quad b.quad <- -matrix(c(12,10),2,1) A.cons <- matrix(c(1,-1,0,1,0,-1),3,2) A.cons b.cons <- c(3.5,0,0) xx <- seq(-1,4,by=0.01) yy <- seq(-1,4,by=0.01) zz <- outer(xx,yy,function(x,y) 2*x^2+x*y+y^2-12*x-10*y) contour(xx,yy,zz,nlev=30) abline(3.5,-1,col=2) abline(h=0,col=2) abline(v=0,col=2) x <- seq(0,3.5,length=20) segments(x,rep(0,length(x)),x,3.5-x,col=2) # initiate: n=0 b.0 <- rep(0,3) x0 <- c(0,0) W0 <- 2:3 # 2. qp <- QP(Q.quad,b.quad+Q.quad%*%x0,A.cons[W0,],b.0[W0]) # 4. qp W1 <- 3 x1 <- x0 points(x1[1],x1[2],col="blue",cex=2) # n=1 # 2. qp <- QP(Q.quad,b.quad+Q.quad%*%x1,A.cons[W1,],b.0[W1]) d1 <- qp$x # 3. alpha <- (b.cons - A.cons%*%x1)/A.cons%*%d1 alpha[A.cons%*%d1 > 0] # 4. qp$lam W2 <- 0 x2 <- x1+d1 points(x2[1],x2[2],col="blue",cex=2) segments(x1[1],x1[2],x2[1],x2[2],lty=2,col="blue",lwd=2) # n=2 # 2. qp <- QP(Q.quad,b.quad+Q.quad%*%x2,A.cons[W2,],b.0[W2]) d2 <- qp$x # 3. alpha <- (b.cons - A.cons%*%x2)/A.cons%*%d2 alpha alpha[A.cons%*%d2 > 0] W3 <- c(1) x3 <- x2+min(alpha[A.cons%*%d2 > 0])*d2 points(x3[1],x3[2],col="blue",cex=2) segments(x2[1],x2[2],x3[1],x3[2],lty=2,col="blue",lwd=2) # n=3 # 2. qp <- QP(Q.quad,b.quad+Q.quad%*%x3,A.cons[W3,],b.0[W3]) d3 <- qp$x # 3. alpha <- (b.cons - A.cons%*%x3)/A.cons%*%d3 alpha alpha[A.cons%*%d3 > 0] A.cons%*%d3 # 4. qp$lam x4 <- x3 + d3 points(x4[1],x4[2],col="blue",cex=2) segments(x3[1],x3[2],x4[1],x4[2],lty=2,col="blue",lwd=2)