# How do Lawson and Hanson solve the unconstrained least squares problem?

For the non-negative least squares problem min(b - a %*% x) subject to x >= 0, the Lawson-Hanson fortran77 implementation (fortran code, R package) often gives a different solution than modern methods for the same set of unconstrained variables.

Modern benchmark methods might include multiway::fnnls, sequential coordinate descent nnls from NNLM, and TNT-NN. Note that bounded Variable Least Squares (bvls::bvls) uses Lawson/Hanson source code and often returns the same results.

Technical background: Non-negative least squares problems are comprised of two independent subproblems:

1. determining what variables will be constrained to zero (the active set)
2. solving the least squares solution for all unconstrained variables

This question deals with the second subproblem. However, note that these subproblems are blurred in coordinate descent which manipulates b to minimize error rather than active set selection. However, NNLM::nnlm is unable to replicate LH77 results and uses gradient descent.

Functions

library(nnls)
library(multiway)

# my own take on tnt-nn
fastnnls <- function(a, b){
x <- solve(a, b)
while(any(x < 0)){
ind <- which(x > 0)
x <- rep(0, length(b))
x[ind] <- solve(a[ind, ind], b[ind])
}
as.vector(x)
}


Reproducible setup

set.seed(123)
X <- matrix(rnorm(2000),100,20)
y <- X %*% runif(20) + rnorm(100)*5
a <- crossprod(X)
b <- as.vector(crossprod(X, y))


Calculate solutions and residuals

x_fnnls <- as.vector(multiway::fnnls(a, b))
x_lh <- nnls::nnls(a, b)$x x_fastnnls <- fastnnls(a, b) all.equal(x_fnnls, x_fastnnls) # [1] TRUE all.equal(x_fnnls, x_lh) # [1] "Mean relative difference: 0.08109434" all.equal(x_fastnnls, x_lh) # [1] "Mean relative difference: 0.08109434"  Note that multiway::fnnls is theoretically supposed to give the same result as the Lawson-Hanson implementation, but rarely does so. However, it does find the same active set, but then the solution is different -- and often worse. Residuals # using mean l2-norm mean(abs(b - as.vector(a %*% x_fnnls))) # [1] 6.515709 mean(abs(b - as.vector(a %*% x_fastnnls))) # [1] 6.515709 mean(abs(b - as.vector(a %*% x_lh))) # [1] 9.086841  In the case of this simulated random example, the error of the Lawson-Hanson algorithm is greater, but I've often seen it significantly lower than other methods for the same active set. So, why is this? What about the LH77 solver gives different results from a modern Cholesky symmetric positive definite solver that is used in R base::solve, Armadillo arma::solve, etc.? I have tried to read the Fortran code and had to throw in the towel. Thoughts: • Does this come down to the method they use to solve the equations? For instance, do they use QR instead of Cholesky decomposition -> forward -> backward substitution, and would this matter? • Do LH77 use gradient optimization after finding the final active set, and if so, why are the coordinate descent solvers in NNLM::nnlm unable to match the performance of LH77 in well-determined systems? • It appears that sometimes the LH77 active set is actually different as well, and this leads to a slightly different (and almost invariably better) result. Thus, this difference seems to be a feature of their solver entirely. • Is there some NNLS-specific method for solving least squares equations that Lawson/Hanson discovered that all follow-up studies failed to replicate, because we like out-of-the-box solvers? • Is the true solution of NNLS actually not solvable by coordinate descent, and the best we can get is an approximation while LH77 somehow finds that true minimum? This is a loaded question, I know, but I'm betting my reputation on it that somebody somewhere around here has a rabbit in their hat. ## 1 Answer ### 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 vector $$Y$$, the other function requires the matrix $$X^TX$$ and vector $$X^Ty$$. 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 } else { new_active <- new_s == 0 ### in the for loop we had only been decreasing the active set ### but the feasible_set function also increases the active it and we need to adapt accrdingly } 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

• Impressive effort, thank you for your time! I wasn't taking XtX and Xty, that's all! Now fastnnls and multiway are giving the same results as nnls::nnls. Just for the record, the fastnnls algorithm that I provided above is only the first part of my method :) The second stage is coordinate descent which achieves the same result as TNT-NN without manipulating the active set. My understanding of TNT-NN is also that they try to manipulate the active set after each iteration, which ultimately slows it down a lot. Mar 2, 2021 at 15:10
• And I would be remiss to not acknowledge your excellent and interpretable implementation of the TNT-NN method, even if it's not the most efficient approach to the problem. Mar 2, 2021 at 15:28
• @zdebruine thank you for your comments about the interpretability. It made me look again at my code and I found a tiny error which I have corrected now. I can not guarantee whether there aren't more errors. Mar 2, 2021 at 16:04