# R/auxillary.R In isotone: Active Set and Generalized PAVA for Isotone Optimization

#### Defines functions is.negis.posis.activeweightedMidRangeweightedFractileweightedMediancheckSolmkLagrangemkIndiwarshallb2ataTxxTaTx

```#This file contains various auxillary functions

#reorganize constraints matrix according to the predictors and take the difference
#e.g. in a total order the difference is 1
aTx <- function(a, x) {
if (is.vector(x)) return(x[a[,1]]-x[a[,2]])
return(x[a[,1],]-x[a[,2],])
}

xT<-function(x) {
if (is.vector(x)) return(as.matrix(x))
else return(t(x))
}

#computes A'lambda which should be equal to gradient
taTx<-function(a,x,n) {
m <- nrow(a)
h <- rep(0,n)
for (i in 1:m) {
h[a[i,1]] <- h[a[i,1]]+x[i]
h[a[i,2]] <- h[a[i,2]]-x[i]
}
return(h)
}

#generates matrix A(I) (i.e. constraint matrix of active sets) as defined in the paper (with 1, -1, and 0)
b2a <- function(b,n) {
m <- nrow(b)                 #number of active sets
a <- matrix(0,m,n)
for (i in 1:m) {             #runs over active sets
a[i,b[i,1]]<- 1            #1 and -1 depending on constraint specification
a[i,b[i,2]]<- -1
}
return(a)
}

#Warshall algorithm for transitive closure
warshall <- function(a) {
n <- nrow(a)
for (j in 1:n) {
for(i in 1:n) {
if (a[i,j]==1) a[i,]<-pmax(a[i,],a[j,])
}
}
return(a)
}

mkIndi<-function(a, n) {
im <- matrix(0, n, n)
m <- nrow(a)
for (i in 1:m) im[a[i,1],a[i,2]] <- im[a[i,2],a[i,1]] <- 1
im <- im + diag(n)
return(t(unique(warshall(im))))
}

#Computes the KKT vector (Lagrange multiplier lambda)
mkLagrange<-function(b, g) {
ta <- t(b2a(b, length(g)))           #matrix A(I) (1, -1, 0) for active set constraints
qa <- qr(ta)                         #QR decomposition
return(qr.coef(qa,g))                #returns KKT vector lambda (solves A*lambda=gradient; stationarity condition)
}

#checks KKT conditions for isotonicity
checkSol<-function(y, gy, a, ay, hl, lbd, ups) {
ckFeasibility <- min(ay)            #Ax >= 0 (primal feasibility)
ckLagrange <- min(lbd)              #lambda >= 0 (dual feasibility)
ckCompSlack <- sum(ay*lbd)          #lambda'Ax = 0 (complementary slackness)
}

weightedMedian<-function(x,w=rep(1,length(x))){
ox<-order(x); x<-x[ox]; w<-w[ox]; k<-1
low<-cumsum(c(0,w)); up<-sum(w)-low; df<-low-up
repeat{
if (df[k] < 0) k<-k+1
else if (df[k] == 0) return((w[k]*x[k]+w[k-1]*x[k-1])/(w[k]+w[k-1]))
else return(x[k-1])
}
}

weightedFractile<-function(x,w=rep(1,length(x)),a=1,b=1){
ox<-order(x); x<-x[ox]; w<-w[ox]; k<-1
low<-cumsum(c(0,w)); up<-sum(w)-low; df<-a*low-b*up
repeat{
if (df[k] < 0) k<-k+1
else if (df[k] == 0) return((w[k]*x[k]+w[k-1]*x[k-1])/(w[k]+w[k-1]))
else return(x[k-1])
}
}

weightedMidRange<-function(x,w=rep(1,length(x))){
s<-0; n<-length(x)
if (n==1) return(x)
for (i in 1:(n-1)) for(j in (i+1):n) {
t<-w[i]*w[j]*abs(x[i]-x[j])/(w[i]+w[j])
if (t > s) {
s<-t; i0<-i; j0<-j
}
}
return((w[i0]*x[i0]+w[j0]*x[j0])/(w[i0]+w[j0]))
}

#which constraint is active (inactives are 0)
is.active <- function(f,ups=1e-12) which(abs(f) < ups)

is.pos<-function(x,ups=1e-12) x > -ups

is.neg<-function(x,ups=1e-12) x < ups
```

## Try the isotone package in your browser

Any scripts or data that you put into this service are public.

isotone documentation built on May 2, 2019, 4:41 p.m.