swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: fdf16693b000f3e309c56091892c61f0ec9fd670 authored by Hans W. Borchers on 21 November 2017, 16:15:00 UTC
version 2.1.1
version 2.1.1
Tip revision: fdf1669
lsqnonneg.R
lsqnonneg <- function(C, d) {
stopifnot(is.numeric(C), is.numeric(d))
if (!is.matrix(C) || !is.vector(d))
stop("Argument 'C' must be a matrix, 'd' a vector.")
m <- nrow(C); n <- ncol(C)
if (m != length(d))
stop("Arguments 'C' and 'd' have nonconformable dimensions.")
tol = 10 * eps() * norm(C, type = "2") * (max(n, m) + 1)
x <- rep(0, n) # initial point
P <- logical(n); Z <- !P # non-active / active columns
resid <- d - C %*% x
w <- t(C) %*% resid
wz <- numeric(n)
# iteration parameters
outeriter <- 0; it <- 0
itmax <- 3 * n; exitflag <- 1
while (any(Z) && any(w[Z] > tol)) {
outeriter <- outeriter + 1
z <- numeric(n)
wz <- rep(-Inf, n)
wz[Z] <- w[Z]
im <- which.max(wz)
P[im] <- TRUE; Z[im] <- FALSE
z[P] <- qr.solve(C[, P], d)
while (any(z[P] <= 0)) {
it <- it + 1
if (it > itmax) stop("Iteration count exceeded.")
Q <- (z <= 0) & P
alpha <- min(x[Q] / (x[Q] - z[Q]))
x <- x + alpha*(z - x)
Z <- ((abs(x) < tol) & P) | Z
P <- !Z
z <- numeric(n)
z[P] <- qr.solve(C[, P], d)
}
x <- z
resid <- d - C %*% x
w <- t(C) %*% resid
}
return(list(x = x, resid.norm = sum(resid*resid)))
}