Nothing
make.ids <- structure(function
### Construct matrices and vectors with their unique ids
(...
### Named lists of integer vectors of length <= 2. Each item
### represents an optimization variable and its dimension.
){
vars <- list(...)
ids <- list()
i <- 1
for(N in names(vars)){
d <- vars[[N]]
n <- prod(d)
end <- i+n-1
a <- structure(array(i:end,d),class="id",variable=N)
i <- end+1
ids[[N]] <- a
}
ids
### Named list of optimization variable ids, which will be used to
### construct the standard form matrices.
},ex=function(){
vars <- make.ids(xi=5,beta0=1,beta=2)
with(vars, xi >= 0)
})
">=.id" <- function(a,constant){
### makes a list of inequality constraints, one for each variable
a[] >= constant
}
"*.id" <- function(a,constant){
### multiply coefficients, returning a linear combination.
sum(a[]) * constant
}
"*.lincomb" <- function(vars,x){
### multiply coefficients, returning a linear combination
if(length(x)==1)x <- rep(x,length(vars))
stopifnot(length(vars)==length(x))
for(i in seq_along(vars)){
vars[[i]]$coef <- vars[[i]]$coef * x[i]
}
structure(vars,class="lincomb")
}
"*.lincomblist" <- function(vars,x){
### multiply coefficients, returning a linear combination list
stopifnot(is.numeric(x))
if(length(x)==1)x <- rep(x,length(vars))
stopifnot(length(x)==length(vars))
for(i in seq_along(vars)){
vars[[i]] <- vars[[i]] * x[i]
}
vars
}
sum.id <- function(a,...){
### sum all variables, returning a linear combination
sum(a[])
}
">=.lincomb" <- function(vars,constant){
### make an inequality constraint from a linear combination
arglist <- c(list("geq",constant),vars)
do.call(constrain,arglist)
}
"==.lincomb" <- function(vars,constant){
### make an equality constraint from a linear combination
arglist <- c(list("eq",constant),vars)
do.call(constrain,arglist)
}
">=.lincomblist" <- function(lcs,constant){
### make a list of inequality constraints from a list of linear combinations
L <- lapply(lcs,function(lc)lc>=constant)
if(length(L)==1)L <- L[[1]]
L
}
sum.lincomblist <- function(lcs,...){
### collapse a list of linear combinations into a single linear combination
structure(lapply(lcs,"[[",1),class="lincomb")
}
"+.lincomb" <- function(v1,v2){
### combine optimization variables in a linear combination
if(class(v2)%in%c("id","lincomblist"))v2 <- sum(v2)
structure(c(v1,v2),class="lincomb")
}
### Model: define a less than operator for a linear combination of
### optimization variables, which returns a constraint. Define an
### addition operator for a linear combination, which returns another
### linear combination, which is just a list of variables and
### coefficients.
"[.id" <- function(var,i,j){
### Make some variables from an id matrix.
N <- attr(var,"variable")
linear.combinations <- list()
if(missing(i))i <- 1:nrow(var)
if(length(dim(var))==1){
for(ii in i){
lc <- structure(list(variable(N,ii)),class="lincomb")
linear.combinations <- c(linear.combinations,list(lc))
}
}else{
if(missing(j))j <- 1:ncol(var)
for(ii in i)for(jj in j){
lc <- structure(list(variable(N,ii,jj)),class="lincomb")
linear.combinations <- c(linear.combinations,list(lc))
}
}
structure(linear.combinations,class="lincomblist")
}
constrain <- structure(function
### Make a constraint list. Write each constraint as c_1 v_1 + ... +
### c_n v_n (in)equality constant, where c_i are real coefficients and
### v_i are optimization variables.
(equality,
### "eq" or "geq" designating whether this is an equality or an
### inequality constraint.
constant,
### constant on the right side of the inequality.
...
### Variables to include in this constraint, usually using something
### like variable("xi",i,j,-1).
){
variables <- list(...)
vars <- lapply(variables,structure,class=c("variable","list"))
structure(list(variables=vars,equality=equality,constant=constant),
class=c("constraint","list"))
},ex=function(){
constraints <- list()
for(i in 1:3){
constraints <-
c(constraints,list(constrain("geq",0,variable("xi",i))))
}
lapply(1:3,function(i)constrain("geq",0,variable("xi",i)))
})
variable <- structure(function
### Make a list that represents an optimization variable and its
### coefficient in the context of a constraint.
(name,
### Name of the optimization variable, created with make.ids.
i=1,
### If the optimization variable is a vector, the element. If it is a
### matrix, the row.
j=NULL,
### Column of the optimization variable.
coef=1
### Coefficient of the optimization variable.
){
structure(list(variable=name,i=i,j=j,coef=coef),
class=c("variable","list"))
### List with elements "variable" "i" "j" "coef", of class "variable".
},ex=function(){
variable("xi",2,coef=-1)
variable("alpha",3,3)
})
print.variable <- function(x,...){
### print method that shows variable indices and coefficient
index <- if(is.null(x$j)){
sprintf("%d",x$i)
}else{
sprintf("%d,%d",x$i,x$j)
}
cat(sprintf("%10.3f*%s[%s]\n",x$coef,x$variable,index))
}
print.constraint <- function(x,...){
### print method that shows each variable and constant
for(v in x$variables){
print(v)
}
key <- c(geq=">=",eq="==")
cat(key[x$equality],x$constant,"\n")
}
is.constraint <- function
### Check if an object is a constraint.
(x
### object to check.
){
"constraint" %in% class(x)
### TRUE if x is of class "constraint".
}
### Reality checks for constraints.
validate.constraint <- function(x,vars){
problem <- function(msg){
str(x)
cat("invalid constraint: ",msg,".\n",sep="")
stop("constraints did not validate.")
}
if(!is.constraint(x))problem('not classed "constraint"')
if( !(x$equality %in% c("geq","eq")) ){
problem('equality should be "eq" or "geq"')
}
if(!is.numeric(x$constant))problem('constant must be numeric')
if(length(x$constant) != 1)problem('constant must be length 1')
if(!is.list(x$variables))problem('variables must be a list')
v.prob <- function(msg){
str(vars)
problem(msg)
}
for(v in x$variables){
if(!is.character(v$variable))v.prob('variable name not character')
if(length(v$variable) != 1)v.prob('variable name must be length 1')
if(! v$variable %in% names(vars))v.prob('variable name not in var list')
## TODO: verify indices, etc.
}
}
standard.form.constraints <- structure(function
### Convert constraints to standard form matrix and vector.
(constraints,
### List of constraints created using constrain().
ids
### List of optimization variable ids created using make.ids().
){
## check if all constraints are valid wrt ids.
for(constraint in constraints){
validate.constraint(constraint,ids)
}
## sort by type of constraint:
eq <- sapply(constraints,"[[","equality")
meq <- sum(eq=="eq") ## number of equality constraints
ordered.constraints <- constraints[order(eq)]
b0 <- sapply(ordered.constraints,"[[","constant")
n.constraints <- length(ordered.constraints)
n.variables <- length(unlist(ids))
A <- matrix(0,n.variables,n.constraints)
for(cn in seq_along(ordered.constraints)){
constraint <- ordered.constraints[[cn]]
for(v in constraint$variables){
a <- unclass(ids[[v$variable]])
if(is.null(a))stop(sprintf("variable %s not found",v$variable))
A.row <- if(is.null(v$j)){
a[v$i]
}else{
a[v$i,v$j]
}
A[A.row,cn] <- v$coef
}
}
list(A=A,b0=b0,meq=meq)
### List of variable coefficients and constraint constants.
},ex=function(){
## example: optimization variable \alpha\in\RR^3 subject to the
## constraint that it must be on the standard simplex:
opt.vars <- make.ids(alpha=3)
constraints <- with(opt.vars,c(alpha >= 0,list(sum(alpha) == 1)))
standard.form.constraints(constraints,opt.vars)
## linear svm example
set.seed(1)
p <- 2
y <- rep(c(-1,1),each=20)
x <- replicate(p,rnorm(length(y),y))
plot(x,col=y+2,asp=1)
n <- nrow(x)
vars <- make.ids(slack=n,intercept=1,normal=p)
constraints <- vars$slack >= 0
for(i in 1:n){
ivars <- with(vars,intercept*y[i] + sum(normal)*(x[i,]*y[i]) + slack[i])
constraints <- c(constraints,list(ivars >= 1))
}
solver.args <- standard.form.constraints(constraints,vars)
n.vars <- length(unlist(vars))
Dvec <- rep(1e-6,n.vars)
Dvec[vars$normal] <- 1
D <- diag(Dvec)
d <- rep(0,n.vars)
d[vars$slack] <- -1 ## C == 1
sol <- quadprog::solve.QP(D,d,solver.args$A,solver.args$b0)
slack <- sol$solution[vars$slack]
normal <- sol$solution[vars$normal]
intercept <- sol$solution[vars$intercept]
title(paste("A linear Support Vector Machine (SVM):",
"margin SVs circled, slack drawn in red for other SVs"))
abline(-intercept/normal[2],-normal[1]/normal[2])
abline((1-intercept)/normal[2],-normal[1]/normal[2],lty="dotted")
abline((-1-intercept)/normal[2],-normal[1]/normal[2],lty="dotted")
f <- function(x)intercept+sum(normal*x)
yfx <- apply(x,1,f)*y
on.margin <- abs(yfx-1)<1e-6
points(x[on.margin,],cex=2,col=y[on.margin]+2)
i <- yfx<1
## these 2 complicated formulas calculate the point on margin where
## the data point starts picking up slack
x1 <- ((y[i]-intercept)*normal[1]-
normal[1]*normal[2]*x[i,2]+
normal[2]^2*x[i,1])/(normal[2]^2+normal[1]^2)
x2 <- (y[i]-intercept)/normal[2]-normal[1]/normal[2]*x1
segments(x[i,1],x[i,2],x1,x2,col="red")
l2norm <- function(x)sqrt(sum(x^2))
rbind(apply(cbind(x1,x2)-x[i,],1,l2norm)*l2norm(normal),slack[i])
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.