#--------------------------------CUT HERE------------------------------------- #! /bin/sh # # This is a shell archive. Save this into a file, edit it # and delete all lines above this comment. Then give this # file to sh by executing the command "sh file". The files # will be extracted into the current directory owned by # you with default permissions. # # The files contained herein are: # # -rw-r--r-- 1 colin user 2535 Feb 27 14:15 Fqpns # -rw-r--r-- 1 colin user 2656 Feb 27 14:15 Oqpns # -rw-r--r-- 1 colin user 563 Feb 27 14:15 Qqpns # -rw-r--r-- 1 colin user 202 Feb 27 14:15 README # echo 'x - Fqpns' if test -f Fqpns; then echo 'shar: not overwriting Fqpns'; else sed 's/^X//' << '________This_Is_The_END________' > Fqpns Xqpns <- function(A, b, Q, cc, tol=1e-12){ X# X# Solve nonsingular quadratic programming problem X# X# find x to minimize -x^t Q x/2 + c^t x s.t. Ax + S = b and S >= 0 X# thus solve Qx + A^t a = c; Ax + S = b s.t. a >= 0, S >= 0 and a^t S = 0 X# let QQ = Q^{-1}; substitute x = QQ ( c - A^t a ), X# and thus solve -A QQ A^t a + s = b - A QQ c (with the above constraints) X# X# EXAMPLE 1: choose x to minimize x^t x s.t. Ax >= b: X# set A = A, b = b, Q = 1, and c = 0 X# X# EXAMPLE 2: choose x to minimize sum(x_i - z_i)^2 s.t. Ax >= 0: X# set A = A, b = 0, Q = 1, and c = z X# X# step 0: setup, follow algebra above, then X# compute slacks S with dual variables a all zero X# sel is a logical vector with T for a non-basic slack, and F otherwise X q <- nrow(A) X zeros <- rep(0,q) X# browser() X if(missing(cc)||(length(cc)==1&&cc==0)) cc <- rep(0,ncol(A)) X if(missing(Q)|Q==1){ X QQc <- cc ; QQAt <- t(A) } else { X QQc <- solve(Q,cc); QQAt <- solve(Q,t(A)) } X D <- - A %*% QQAt X z <- b - A %*% QQc X a <- zeros X neg <- rep(F,q) X S <- z X# browser() X while(any(S < -tol) | any(neg)) { X# step 1: non-basic slacks include non-positive slacks and X# basic dual variables with negative values X# cat("step 1: =",sum(a*S),"\n") X# sel <- S < -tol | a > tol X sel <- S < tol & !neg X cat("basic dual variables: ", seq(q)[sel], "\n") X# step 2: set non-basic slacks to zero X a <- S <- zeros X if(sum(sel) > 0) X a[sel] <- solve(D[sel, sel], z[sel]) X# step 3: given dual variables (a) compute slacks (S) X# an alternative formula is: S1 <- z - D %*% a X if(sum(sel)>0 & sum(!sel)>0) X S[!sel] <- -D[!sel,sel] %*% a[sel] + z[!sel] X# cat("step 3: =",sum(a*S),"\n") X# step 4: identify basic dual variables with negative values to make non basic X neg <- a< -tol X if(any(neg)) cat("basic dual variables made non basic: ", X seq(q)[neg], "\n") X# browser() X# step 5: end iterations if no further change and output x X } X drop(QQc - QQAt %*% a) X} X X X"pascal" <- function(length.wanted,x=1) Xif(length(x)>=length.wanted) x else pascal(length.wanted,c(0,x)+c(x,0)) X"diffs" <- function(n){ X d <- pascal(n+1) X d[seq(2,n+1,2)] <- -d[seq(2,n+1,2)] X d X } X"offdiag" <- function(vec,nc,nr=nc-length(vec)+1){ X# put elements of vec on rows of nr x nc matrix, starting at the diagonals X if(length(nc)>1) nc <- length(nc) X mat <- matrix(0,nr,nc) X rowno <- row(mat) X colno <- col(mat) X for(i in seq(length(vec))) mat[colno-rowno==i-1] <- vec[i] X mat X } X"Asmooth" <- function(z) rbind(offdiag(diffs(1),z), offdiag(diffs(3),z)) X ________This_Is_The_END________ if test `wc -l < Fqpns` -ne 76; then echo 'shar: Fqpns was damaged during transit (should have been 76 bytes)' fi fi ; : end of overwriting check echo 'x - Oqpns' if test -f Oqpns; then echo 'shar: not overwriting Oqpns'; else sed 's/^X//' << '________This_Is_The_END________' > Oqpns X> > # Example 3.1 X> Q31 <- 1 X> c31 <- 0 X> A31 <- matrix(c(-1, -1, -1, -1, -2, 1), 2, 3, byrow = T) X> b31 <- c(-1, -2) X> qpns(A31, b31, Q31, c31) Xbasic dual variables: 1 2 X 0.4285714 0.7142857 -0.1428571 X> > # Example 4.1 X> z <- c(0, 3, 4, 6.5, 8, 10.5, 12) X> pascal(0) X 1 X> pascal(2) X 1 1 X> pascal(9) X 1 8 28 56 70 56 28 8 1 X> diffs(1) X 1 -1 X> diffs(3) X 1 -3 3 -1 X> offdiag(diffs(3),z) X [,1] [,2] [,3] [,4] [,5] [,6] [,7] X[1,] 1 -3 3 -1 0 0 0 X[2,] 0 1 -3 3 -1 0 0 X[3,] 0 0 1 -3 3 -1 0 X[4,] 0 0 0 1 -3 3 -1 X> A <- Asmooth(z) X> A X [,1] [,2] [,3] [,4] [,5] [,6] [,7] X [1,] 1 -1 0 0 0 0 0 X [2,] 0 1 -1 0 0 0 0 X [3,] 0 0 1 -1 0 0 0 X [4,] 0 0 0 1 -1 0 0 X [5,] 0 0 0 0 1 -1 0 X [6,] 0 0 0 0 0 1 -1 X [7,] 1 -3 3 -1 0 0 0 X [8,] 0 1 -3 3 -1 0 0 X [9,] 0 0 1 -3 3 -1 0 X[10,] 0 0 0 1 -3 3 -1 X> z X 0.0 3.0 4.0 6.5 8.0 10.5 12.0 X> qpns(A,0,1,z) Xbasic dual variables: 8 10 Xbasic dual variables: 8 9 10 X 0.000000 2.767857 4.475000 6.271429 8.157143 10.132143 12.196429 X> > # Other smoothing examples X> zz <- seq(10) X> AA <- Asmooth(zz) X> qpns(AA,0,1,zz) X  1 2 3 4 5 6 7 8 9 10 X> zzz <- seq(10) + rnorm(10) X> zzz X  1.009556 1.763031 2.950119 4.919003 5.539653 4.099332 6.451911 6.682013 X  8.389683 9.388330 X> qpns(AA,0,1,zzz) Xbasic dual variables: 5 11 12 14 16 Xbasic dual variables: 5 10 11 12 14 15 16 Xbasic dual variables made non basic: 5 Xbasic dual variables: 10 11 12 14 15 16 X  0.680858 2.186755 3.398901 4.317296 4.941940 5.272833 5.903416 6.833687 X  8.063648 9.593298 X> zzzz <- rnorm(10) X> zzzz X  0.71716396 0.83544225 -0.01712495 1.11915948 0.83364126 -0.86914788 X  1.02496840 0.59624871 -0.14620412 -0.75779374 X> qpns(AA,0,1,zzzz) # there is unlikely to be a solution to this Xbasic dual variables: 2 4 5 7 8 9 11 14 Xbasic dual variables made non basic: 11 14 Xbasic dual variables: 1 2 3 4 5 6 7 8 9 10 12 13 15 16 XError in solve.qr(a, b): apparently singular matrix XDumped X> qpns(AA,0,1,sort(zzzz)) Xbasic dual variables: 10 12 14 16 Xbasic dual variables made non basic: 14 Xbasic dual variables: 10 11 12 13 16 Xbasic dual variables: 7 10 11 12 13 14 16 Xbasic dual variables made non basic: 7 Xbasic dual variables: 10 11 12 13 14 16 X  -1.0111084 -0.5668354 -0.1789004 0.1526964 0.4279552 0.6468759 X  0.8094584 0.9157029 1.0195398 1.1209690 X> > > ________This_Is_The_END________ if test `wc -l < Oqpns` -ne 81; then echo 'shar: Oqpns was damaged during transit (should have been 81 bytes)' fi fi ; : end of overwriting check echo 'x - Qqpns' if test -f Qqpns; then echo 'shar: not overwriting Qqpns'; else sed 's/^X//' << '________This_Is_The_END________' > Qqpns Xsource("Fqpns") Xoptions(echo=T) X X# Example 3.1 XQ31 <- 1 Xc31 <- 0 XA31 <- matrix(c(-1, -1, -1, -1, -2, 1), 2, 3, byrow = T) Xb31 <- c(-1, -2) Xqpns(A31, b31, Q31, c31) X X# Example 4.1 Xz <- c(0, 3, 4, 6.5, 8, 10.5, 12) Xpascal(0) Xpascal(2) Xpascal(9) Xdiffs(1) Xdiffs(3) Xoffdiag(diffs(3),z) XA <- Asmooth(z) XA Xz Xqpns(A,0,1,z) X X# Other smoothing examples Xzz <- seq(10) XAA <- Asmooth(zz) Xqpns(AA,0,1,zz) Xzzz <- seq(10) + rnorm(10) Xzzz Xqpns(AA,0,1,zzz) Xzzzz <- rnorm(10) Xzzzz Xqpns(AA,0,1,zzzz) # there is unlikely to be a solution to this Xqpns(AA,0,1,sort(zzzz)) X X ________This_Is_The_END________ if test `wc -l < Qqpns` -ne 36; then echo 'shar: Qqpns was damaged during transit (should have been 36 bytes)' fi fi ; : end of overwriting check echo 'x - README' if test -f README; then echo 'shar: not overwriting README'; else sed 's/^X//' << '________This_Is_The_END________' > README XS functions for quadratic programming: see Goodall and Brown (1994), XIn: A Probability, Statistics and Optimization, A Tribute Xto Peter Whittle (Wiley). Contributed by Colin Goodall (crg@stat.psu.edu). ________This_Is_The_END________ if test `wc -l < README` -ne 3; then echo 'shar: README was damaged during transit (should have been 3 bytes)' fi fi ; : end of overwriting check exit 0