# R/linprogPD.R In clime: Constrained L1-Minimization for Inverse (Covariance) Matrix Estimation

#### Defines functions linprogPD

linprogPD <- function(x0, A, b, epsilon, pdtol=1e-3, pdmaxiter=50) {
## Solves
## min ||x||_1   subject to  ||Ax-b||_\infty <= epsilon
## Adapted from Matlab code for Dantzig Selector by J. Romberg
N <- length(x0)
x0 <- matrix(x0, nrow=N)
b <- matrix(b, nrow=N)

alpha=0.01;
beta=0.5;
mu = 10;

if (max(abs(A%*%x0 - b)) > epsilon) {
stop("Infeasible starting point!")
}

x <- x0
u <- 0.95*abs(x0) + 0.1*max(abs(x0))

Atr <- A%*%x - b

fu1 <- x - u
fu2 <- -x - u
fe1 <- Atr - epsilon
fe2 <- -Atr - epsilon
lamu1 <- -1/fu1
lamu2 <- -1/fu2
lame1 <- -1/fe1
lame2 <- -1/fe2

AtAv <- t(A)%*%(lame1 - lame2)

sdg <- - sum(c(fu1, fu2, fe1, fe2)* c(lamu1, lamu2, lame1, lame2))

tau <- mu*(4*N)/sdg

rdual <- gradf0 + c(lamu1 - lamu2 + AtAv, -lamu1 - lamu2)
rcent <- -c(lamu1*fu1, lamu2*fu2, lame1*fe1, lame2*fe2) - 1/tau
resnorm <- sqrt(sum(rdual^2,rcent^2))

pditer <- 0
done <- (sdg < pdtol) | (pditer >= pdmaxiter)

while(!done) {
w2 <- -1 - (1/fu1 + 1/fu2)/tau

sig11 <- -lamu1/fu1 - lamu2/fu2
sig12 <- lamu1/fu1 - lamu2/fu2
siga <- -(lame1/fe1 + lame2/fe2)
sigx <- sig11 - sig12^2/sig11

w1 <- -( t(A)%*%(1/fe2 - 1/fe1) + 1/fu2 - 1/fu1)/tau
w1p <- w1 - (sig12/sig11)*w2

Hp <- t(A)%*%diag(as.vector(siga))%*%A + diag(as.vector(sigx))

dx <- solve(Hp, w1p)

if (rcond(Hp) < 1e-14) {
warning("Ill conditioned matrix.  Previous iterate matrix returned! (May increase perturb/lambda.)")
xp <- x
return(xp)
}

du <- w2/sig11 - (sig12/sig11)*dx

dlamu1 <- -(lamu1/fu1)*(dx-du) - lamu1 - 1/(fu1*tau)
dlamu2 <- -(lamu2/fu2)*(-dx-du) - lamu2 - 1/(fu2*tau)

dlame1 <- -(lame1/fe1)*(AtAdx) - lame1 - 1/(fe1*tau)
dlame2 <- (lame2/fe2)*(AtAdx) - lame2 - 1/(fe2*tau)

iu1 <- dlamu1 < 0
iu2 <- dlamu2 < 0
ie1 <- dlame1 < 0
ie2 <- dlame2 < 0
ifu1 <- (dx-du) > 0
ifu2 <- (-dx-du) > 0

smax <- min( -lamu1[iu1]/dlamu1[iu1], -lamu2[iu2]/dlamu2[iu2], -lame1[ie1]/dlame1[ie1], -lame2[ie2]/dlame2[ie2], -fu1[ifu1]/(dx[ifu1] - du[ifu1]), -fu2[ifu2]/(-dx[ifu2] -du[ifu2]), -fe1[ife1]/AtAdx[ife1], -fe2[ife2]/( - AtAdx[ife2])   )
smax <- min(1,smax)

s <- 0.99*smax

suffdec <- FALSE
backiter <- 0

while(!suffdec) {
xp <- x + s*dx
up <- u + s*du

fu1p <- fu1 + s*(dx - du)
fu2p <- fu2 + s*(-dx-du)

lamu1p <- lamu1 + s*dlamu1
lamu2p <- lamu2 + s*dlamu2
lame1p <- lame1 + s*dlame1
lame2p <- lame2 + s*dlame2

rdp <- gradf0 + c(lamu1p - lamu2p + AtAvp, -lamu1p - lamu2p)
rcp <- -c(lamu1p*fu1p, lamu2p*fu2p, lame1p*fe1p, lame2p*fe2p) - 1/tau
suffdec <- sqrt( sum(rdp^2, rcp^2)) < (1 - alpha*s)*resnorm

s <- beta*s
backiter <- backiter+1

if (backiter > 32) {
warning("Backtracking stuck.  Previous iterate matrix returned!")
xp <- x
return(xp)
}
}
x <- xp
u <- up
Atr <- Atrp
AtAv <- AtAvp
fu1 <- fu1p
fu2 <- fu2p
fe1 <- fe1p
fe2 <- fe2p
lamu1 <- lamu1p
lamu2 <- lamu2p
lame1 <- lame1p
lame2 <- lame2p

sdg <- -sum( c(fu1, fu2, fe1, fe2)*c(lamu1, lamu2, lame1, lame2))
tau <- mu*(4*N)/sdg
rdual <- rdp
rcent <- -c(fu1, fu2, fe1, fe2)*c(lamu1, lamu2, lame1, lame2) - 1/tau
resnorm <- sqrt(sum(rdual^2, rcent^2))

pditer <- pditer+1
done <- (sdg < pdtol) | (pditer >= pdmaxiter)
}

return(xp)
}


## Try the clime package in your browser

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

clime documentation built on June 22, 2022, 5:07 p.m.