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

Question

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.

Answer 1

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