c======================================================================= subroutine bvls(key, m, n, a, b, bl, bu, x, w, act, zz, istate, + loopA) c======================================================================= implicit double precision (a-h, o-z) c c$$$$ calls qr c--------------------Bounded Variable Least Squares--------------------- c c Robert L. Parker and Philip B. Stark Version 3/19/90 c c Solves the problem: c c min || a.x - b || such that bl <= x <= bu c 2 c where c x is an unknown n-vector c a is a given m by n matrix c b is a given m-vector c bl is a given n-vector of lower bounds on the c components of x. c bu is a given n-vector of upper bounds on the c components of x. c c c----------------------------------------------------------------------- c Input parameters: c c m, n, a, b, bl, bu see above. Let mm=min(m,n). c c If key = 0, the subroutine solves the problem from scratch. c c If key > 0 the routine initializes using the user's guess about c which components of x are `active', i.e. are stricly within their c bounds, which are at their lower bounds, and which are at their c upper bounds. This information is supplied through the array c istate. istate(n+1) should contain the total number of components c at their bounds (the `bound variables'). The absolute values of the c first nbound=istate(n+1) entries of istate are the indices c of these `bound' components of x. The sign of istate(j), j=1,..., c nbound, indicates whether x(|istate(j)|) is at its upper or lower c bound. istate(j) is positive if the component is at its upper c bound, negative if the component is at its lower bound. c istate(j), j=nbound+1,...,n contain the indices of the components c of x that are active (i.e. are expected to lie strictly within c their bounds). When key > 0, the routine initially sets the active c components to the averages of their upper and lower bounds: c x(j)=(bl(j)+bu(j))/2, for j in the active set. c c----------------------------------------------------------------------- c Output parameters: c c x the solution vector. c c w(1) the minimum 2-norm || a.x-b ||. c c istate vector indicating which components of x are active and c which are at their bounds (see the previous paragraph). c istate can be supplied to the routine to give it a good c starting guess for the solution. c c loopA number of iterations taken in the main loop, Loop A. c c----------------------------------------------------------------------- c Working arrays: c c w dimension n. act dimension m*(mm+2). c zz dimension m. istate dimension n+1. c c----------------------------------------------------------------------- c Method: active variable method along the general plan of NNLS by c Lawson & Hanson, "Solving Least Squares Problems," 1974. See c Algorithm 23.10. Step numbers in comment statements refer to their c scheme. c c----------------------------------------------------------------------- c A number of measures are taken to enhance numerical reliability: c c 1. As noted by Lawson and Hanson, roundoff errors in the computation c of the gradient of the misfit may cause a component on the bounds c to appear to want to become active, yet when the component is added c to the active set, it moves away from the feasible region. In this c case the component is not made active, the gradient of the misfit c with respect to a change in that component is set to zero, and the c program returns to the Kuhn-Tucker test. Flag ifrom5 is used in c this test, which occurs at the end of Step 6. c c c 2. When the least-squares minimizer after Step 6 is infeasible, it c is used in a convex interpolation with the previous solution to c obtain a feasible vector. The constant in this interpolation is c supposed to put at least one component of x on a bound. There can c be difficulties: c c 2a. Sometimes, due to roundoff, no interpolated component ends up on c a bound. The code in Step 11 uses the flag jj, computed in Step 8, c to ensure that at least the component that determined the c interpolation constant alpha is moved to the appropriate bound. c This guarantees that what Lawson and Hanson call `Loop B' is finite. c c 2b. The code in Step 11 also incorporates Lawson and Hanson's feature c that any components remaining infeasible at this stage (which must c be due to roundoff) are moved to their nearer bound. c c c 3. If the columns of a passed to qr are linearly dependent, the new c potentially active component is not introduced: the gradient of the c misfit with respect to that component is set to zero, and control c returns to the Kuhn-Tucker test. c c c 4. When some of the columns of a are approximately linearly c dependent, we have observed cycling of active components: a c component just moved to a bound desires immediately to become c active again; qr allows it to become active and a different c component is moved to its bound. This component immediately wants c to become active, which qr allows, and the original component is c moved back to its bound. We have taken two steps to avoid this c problem: c c 4a. First, the column of the matrix a corresponding to the new c potentially active component is passed to qr as the last column of c its matrix. This ordering tends to make a component recently moved c to a bound fail the test mentioned in (1), above. c c 4b. Second, we have incorporated a test that prohibits short cycles. c If the most recent successful change to the active set was to move c the component x(jj) to a bound, x(jj) is not permitted to reenter c the solution at this stage. This test occurs just after checking c the Kuhn-Tucker conditions, and uses the flag jj, set in Step 8. c The flag jj is reset after Step 6 if Step 6 was entered from c Step 5 indicating that a new component has successfully entered the c active set. The test for resetting jj uses the flag ifrom5, c which will not equal zero in case Step 6 was entered from Step 5. c c dimension a(m,n), b(m), x(n), bl(n), bu(n) c dimension w(n), act(m,min(m,n)+2), zz(m), istate(n+1) dimension w(n), act(m,m+2), zz(m), istate(n+1) c data eps/1.0d-11/ c c----------------------First Executable Statement----------------------- c c Step 1. Initialize everything--active and bound sets, initial c values, etc. c c Initialize flags, etc. mm=min(m,n) mm1 = mm + 1 jj = 0 ifrom5 = 0 c Check consistency of given bounds bl, bu. bdiff = 0.0 do 1005 j=1, n bdiff=max(bdiff, bu(j)-bl(j)) if (bl(j) .gt. bu(j)) then print *,' Inconsistent bounds in BVLS. ' stop endif 1005 continue if (bdiff .eq. 0.0) then print *,' No free variables in BVLS--check input bounds.' stop endif c c In a fresh initialization (key = 0) bind all variables at their lower c bounds. If (key != 0), use the supplied istate vector to c initialize the variables. istate(n+1) contains the number of c bound variables. The absolute values of the first c nbound=istate(n+1) entries of istate are the indices of the bound c variables. The sign of each entry determines whether the indicated c variable is at its upper (positive) or lower (negative) bound. if (key .eq. 0) then nbound=n nact=0 do 1010 j=1, nbound istate(j)=-j 1010 continue else nbound=istate(n+1) endif nact=n - nbound if ( nact .gt. mm ) then print *, ' Too many active variables in BVLS starting solution!' stop endif do 1100 k=1, nbound j=abs(istate(k)) if (istate(k) .lt. 0) x(j)=bl(j) if (istate(k) .gt. 0) x(j)=bu(j) 1100 continue c c In a warm start (key != 0) initialize the active variables to c (bl+bu)/2. This is needed in case the initial qr results in c active variables out-of-bounds and Steps 8-11 get executed the c first time through. do 1150 k=nbound+1,n kk=istate(k) x(kk)=(bu(kk)+bl(kk))/2 1150 continue c c Compute bnorm, the norm of the data vector b, for reference. bsq=0.0 do 1200 i=1, m bsq=bsq + b(i)**2 1200 continue bnorm=sqrt(bsq) c c-----------------------------Main Loop--------------------------------- c c Initialization complete. Begin major loop (Loop A). do 15000 loopA=1, 3*n c c Step 2. c Initialize the negative gradient vector w(*). 2000 obj=0.0 do 2050 j=1, n w(j)=0.0 2050 continue c c Compute the residual vector b-a.x , the negative gradient vector c w(*), and the current objective value obj = || a.x - b ||. c The residual vector is stored in the mm+1'st column of act(*,*). do 2300 i=1, m ri=b(i) do 2100 j=1, n ri=ri - a(i,j)*x(j) 2100 continue obj=obj + ri**2 do 2200 j=1, n w(j)=w(j) + a(i,j)*ri 2200 continue act(i,mm1)=ri 2300 continue c c Converged? Stop if the misfit << || b ||, or if all components are c active (unless this is the first iteration from a warm start). if (sqrt(obj) .le. bnorm*eps .or. + (loopA .gt. 1 .and. nbound .eq. 0)) then istate(n+1)=nbound w(1)=sqrt(obj) return endif c c Add the contribution of the active components back into the residual. do 2500 k=nbound+1, n j=istate(k) do 2400 i=1, m act(i,mm1)=act(i,mm1) + a(i,j)*x(j) 2400 continue 2500 continue c c The first iteration in a warm start requires immediate qr. if (loopA .eq. 1 .and. key .ne. 0) goto 6000 c c Steps 3, 4. c Find the bound element that most wants to be active. 3000 worst=0.0 it=1 do 3100 j=1, nbound ks=abs(istate(j)) bad=w(ks)*sign(1, istate(j)) if (bad .lt. worst) then it=j worst=bad iact=ks endif 3100 continue c c Test whether the Kuhn-Tucker condition is met. if (worst .ge. 0.0 ) then istate(n+1)=nbound w(1)=sqrt(obj) return endif c c The component x(iact) is the one that most wants to become active. c If the last successful change in the active set was to move x(iact) c to a bound, don't let x(iact) in now: set the derivative of the c misfit with respect to x(iact) to zero and return to the Kuhn-Tucker c test. if ( iact .eq. jj ) then w(jj)=0.0 goto 3000 endif c c Step 5. c Undo the effect of the new (potentially) active variable on the c residual vector. if (istate(it) .gt. 0) bound=bu(iact) if (istate(it) .lt. 0) bound=bl(iact) do 5100 i=1, m act(i,mm1)=act(i,mm1) + bound*a(i,iact) 5100 continue c c Set flag ifrom5, indicating that Step 6 was entered from Step 5. c This forms the basis of a test for instability: the gradient c calculation shows that x(iact) wants to join the active set; if c qr puts x(iact) beyond the bound from which it came, the gradient c calculation was in error and the variable should not have been c introduced. ifrom5=istate(it) c c Swap the indices (in istate) of the new active variable and the c rightmost bound variable; `unbind' that location by decrementing c nbound. istate(it)=istate(nbound) nbound=nbound - 1 nact=nact + 1 istate(nbound+1)=iact c if (mm .lt. nact) then print *,' Too many free variables in BVLS!' stop endif c c Step 6. c Load array act with the appropriate columns of a for qr. For c added stability, reverse the column ordering so that the most c recent addition to the active set is in the last column. Also c copy the residual vector from act(., mm1) into act(., mm1+1). 6000 do 6200 i=1, m act(i,mm1+1)=act(i,mm1) do 6100 k=nbound+1, n j=istate(k) act(i,nact+1-k+nbound)=a(i,j) 6100 continue 6200 continue c call qr(m, nact, act, act(1,mm1+1), zz, resq) c c Test for linear dependence in qr, and for an instability that moves c the variable just introduced away from the feasible region c (rather than into the region or all the way through it). c In either case, remove the latest vector introduced from the c active set and adjust the residual vector accordingly. c Set the gradient component (w(iact)) to zero and return to c the Kuhn-Tucker test. if (resq .lt. 0.0 + .or. (ifrom5 .gt. 0 .and. zz(nact) .gt. bu(iact)) + .or. (ifrom5 .lt. 0 .and. zz(nact) .lt. bl(iact))) then nbound=nbound + 1 istate(nbound)=istate(nbound)*sign(1.0d0, x(iact)-bu(iact)) nact=nact - 1 do 6500 i=1, m act(i,mm1)=act(i,mm1) - x(iact)*a(i,iact) 6500 continue ifrom5=0 w(iact)=0.0 goto 3000 endif c c If Step 6 was entered from Step 5 and we are here, a new variable c has been successfully introduced into the active set; the last c variable that was fixed at a bound is again permitted to become c active. if ( ifrom5 .ne. 0 ) jj=0 ifrom5=0 c c Step 7. Check for strict feasibility of the new qr solution. do 7100 k=1, nact k1=k j=istate(k+nbound) if (zz(nact+1-k).lt.bl(j) .or. zz(nact+1-k).gt.bu(j)) goto 8000 7100 continue do 7200 k=1, nact j=istate(k+nbound) x(j)=zz(nact+1-k) 7200 continue c New iterate is feasible; back to the top. goto 15000 c c Steps 8, 9. 8000 alpha=2.0 alf=alpha do 8200 k=k1, nact j=istate(k+nbound) if (zz(nact+1-k) .gt. bu(j)) + alf=(bu(j)-x(j))/(zz(nact+1-k)-x(j)) if (zz(nact+1-k) .lt. bl(j)) + alf=(bl(j)-x(j))/(zz(nact+1-k)-x(j)) if (alf .lt. alpha) then alpha=alf jj=j sj=sign(1.0, zz(nact+1-k)-bl(j)) endif 8200 continue c c Step 10 do 10000 k=1, nact j=istate(k+nbound) x(j)=x(j) + alpha*(zz(nact+1-k)-x(j)) 10000 continue c c Step 11. c Move the variable that determined alpha to the appropriate bound. c (jj is its index; sj is + if zz(jj)> bu(jj), - if zz(jj)