TNT-NN
For the TNT-NN see:
Myre, Joe M., et al. "TNT-NN: a fast active set method for solving large non-negative least squares problems." Procedia Computer Science 108 (2017): 755-764.
https://doi.org/10.1016/j.procs.2017.05.194
To form the active set, TNT-NN first solves an unconstrained least squares problem. Variables that violate the non-negativity constraint are added to the active set. Once a feasible solution is found, where none of the non-negativity constraints are violated, the 2-norm of the residual is used as a measure of fitness and the solution is saved as the current “best” solution.
TNT-NN attempts to modify the active set by iteratively moving some of the variables from the active set back into the unconstrained set. The active set variables are sorted based on their components of the gradient. Variables that show the largest positive gradient components are tested by moving some of them from the active set into the unconstrained set. It is important to note that initially large groups of variables can be moved in a single test. If the new solution
does not improve in fitness, then the solution is rejected and a smaller set of variables is tested. If a group of the variables can be removed from the active set and a new feasible solution is found that is “better”, the solution is saved and the algorithm begins a new iteration. The algorithm reaches convergence when the active set can no longer be modified.
Your implementation is only half part of this solution. It is the part where the new feasible set is searched to test the removal of a variable from the active set. In the code below it is demonstrated. I have copied your fastnnls function and turned it into feasible_set which is only part of the algorithm.
In the article Myre et al speak of " by moving some of them from the active set into the unconstrained set. It is important to note that initially large groups of variables can be moved in a single test". But I am not sure how they do that so in the code, I have been adding them one by one. Probably there are some additional tricks to make faster selections of large groups to be added at once instead of my for loop that tries all variables.
The difference between multiway::fnnls and nnls::nnls
You got a difference because of a small error in your comparison. The one function requires the matrix $X$ and You have used the latter for both functions. In the code below I give an example output.
Example code
### finding the feasible set
feasible_set <- function(a, b, ind){
x <- rep(0, length(b))
x[ind] <- solve(a[ind, ind], b[ind])
while(any(x < 0)){
ind <- which(x > 0)
x <- rep(0, length(b))
x[ind] <- solve(a[ind, ind], b[ind])
}
as.vector(x)
}
### finding the gradients
gradients <- function(b,y,X) {
current_y <- X %*% b
d_y <- y-current_y
gradients <- t(X) %*% d_y
return(gradients)
}
### The algorithm that repeatedly updates the active set
### The updates are done by removing the variable with the highest positive gradient
fastnnls <- function(y,X) {
### Initiation
a <- crossprod(X)
b <- as.vector(crossprod(X, y))
current_active <- rep(TRUE,length(X[1,])) ### start with all variables in active set
current_s <- rep(0,length(X[1,])) ### initial conditions
current_y <- X %*% current_s
current_loss <- sum((y-current_y)^2)
### algorithm that stops untill no improvement can be made
cont <- TRUE
while (cont) {
### add variables based on gradients
### in these four lines the gradients are found and ordered
gradients <- gradients(current_s,y,X)
testing <- which(gradients*current_active>0) ### find out which variables are active and have positive gradients
ord <- order(gradients, decreasing = TRUE)
ord <- ord[ord %in% testing] ### strip the negative or non-active variables
### keep adding variables in a loop while this improves the solution
addition <- 0 ### itterative variable keeping track of the additions
new_active <- current_active
for (i in 1:length(ord)) {
### Try out a new active set with one variable removed
new_active[ord[i]] <- FALSE ### remove 'ord[i]' from active set
new_s <- feasible_set(a,b, ind = which(new_active == FALSE))
new_y <- X %*% new_s
new_loss <- sum((y-new_y)^2)
### Update the solution if the new trial is better
if (new_loss < current_loss) {
addition <- i
current_active <- new_active
current_loss <- new_loss
current_s <- new_s
current_y <- new_y
} else {
break ### skip loop to end
}
}
if (addition == 0) { ### quit while when no addition is made
cont = FALSE
}
if (sum(current_active) == 0) { ### quit if active set is empty (all variables positive)
cont = FALSE
}
### The while loop continues by recomputing the gradients
}
return(current_s)
}
set.seed(123)
X <- matrix(rnorm(2000),100,20)
y <- X %*% runif(20) + rnorm(100)*5
library(nnls)
library(multiway)
data.frame(multiway = multiway::fnnls(a, b),
nnls = nnls::nnls(X, y)$x,
manual = fastnnls(y,X))
Output
> data.frame(multiway = multiway::fnnls(a, b),
+ nnls = nnls::nnls(X, y)$x,
+ manual = fastnnls(y,X))
multiway nnls manual
1 0.610802720 0.610802720 0.610802720
2 0.146121047 0.146121047 0.146121047
3 0.841809005 0.841809005 0.841809005
4 1.131040740 1.131040740 1.131040740
5 0.000000000 0.000000000 0.000000000
6 1.093652478 1.093652478 1.093652478
7 0.725590111 0.725590111 0.725590111
8 0.211525228 0.211525228 0.211525228
9 0.000000000 0.000000000 0.000000000
10 1.472333600 1.472333600 1.472333600
11 0.005740395 0.005740395 0.005740395
12 2.131277775 2.131277775 2.131277775
13 0.000000000 0.000000000 0.000000000
14 0.590923989 0.590923989 0.590923989
15 0.652530944 0.652530944 0.652530944
16 0.717713755 0.717713755 0.717713755
17 1.115162378 1.115162378 1.115162378
18 0.603304661 0.603304661 0.603304661
19 0.000000000 0.000000000 0.000000000
20 0.218073317 0.218073317 0.218073317
