MODULE BoundedLeastSquares
! Summary:
! BVLS solves linear least-squares problems with upper and lower bounds on the
! variables, using an active set strategy. It is documented in the J. of
! Computational Statistics, and can be used iteratively to solve minimum
! l-1, l-2 and l-infinity fitting problems.
! Statement:
! This algorithm may be used freely for non-commercial purposes, and may be
! freely distributed for non-commercial purposes. The authors do not warrant
! the software in any way: use it at your own risk.
! Code converted using TO_F90 by Alan Miller
! Date: 2000-11-23 Time: 23:41:42
IMPLICIT NONE
INTEGER, PARAMETER, PRIVATE :: dp = SELECTED_REAL_KIND(12, 60)
CONTAINS
!=======================================================================
SUBROUTINE bvls(key, m, n, a, b, bl, bu, x, istate, loopa)
!=======================================================================
! N.B. Arguments W, ACT & ZZ have been removed.
INTEGER, INTENT(IN) :: key
INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: a(m,n)
REAL (dp), INTENT(IN OUT) :: b(m)
REAL (dp), INTENT(IN) :: bl(n)
REAL (dp), INTENT(IN) :: bu(n)
REAL (dp), INTENT(OUT) :: x(n)
INTEGER, INTENT(IN OUT) :: istate(n+1)
INTEGER, INTENT(OUT) :: loopa
!--------------------Bounded Variable Least Squares---------------------
! Robert L. Parker and Philip B. Stark Version 3/19/90
! Robert L. Parker Philip B. Stark
! Scripps Institution of Oceanography Department of Statistics
! University of California, San Diego University of California
! La Jolla CA 92093 Berkeley CA 94720-3860
! rlparker@ucsd.edu stark@stat.berkeley.edu
! Copyright of this software is reserved by the authors; however, this
! algorithm and subroutine may be used freely for non-commercial
! purposes, and may be distributed freely for non-commercial purposes.
! The authors do not warrant this software in any way: use it at your
! own risk.
! See the article ``Bounded Variable Least Squares: An Algorithm and
! Applications'' by P.B. Stark and R.L. Parker, in the journal:
! Computational Statistics, vol.10(2), (1995) for further description and
! applications to minimum l-1, l-2 and l-infinity fitting problems, as well
! as finding bounds on linear functionals subject to bounds on variables and
! fitting linear data within l-1, l-2 or l-infinity measures of misfit.
! BVLS solves the problem:
! min || a.x - b || such that bl <= x <= bu
! 2
! where
! x is an unknown n-vector
! a is a given m by n matrix
! b is a given m-vector
! bl is a given n-vector of lower bounds on the components of x.
! bu is a given n-vector of upper bounds on the components of x.
!-----------------------------------------------------------------------
! Input parameters:
! m, n, a, b, bl, bu see above. Let mm=min(m,n).
! If key = 0, the subroutine solves the problem from scratch.
! If key > 0 the routine initializes using the user's guess about which
! components of x are `active', i.e. are stricly within their bounds,
! which are at their lower bounds, and which are at their upper bounds.
! This information is supplied through the array istate. istate(n+1)
! should contain the total number of components at their bounds (the
! `bound variables'). The absolute values of the first nbound=istate(n+1)
! entries of istate are the indices of these `bound' components of x.
! The sign of istate(j), j=1,..., nbound, indicates whether x(|istate(j)|)
! is at its upper or lower bound. istate(j) is positive if the component
! is at its upper bound, negative if the component is at its lower bound.
! istate(j), j=nbound+1,...,n contain the indices of the components of x
! that are active (i.e. are expected to lie strictly within their bounds).
! When key > 0, the routine initially sets the active components to the
! averages of their upper and lower bounds:
! x(j)=(bl(j)+bu(j))/2, for j in the active set.
!-----------------------------------------------------------------------
! Output parameters:
! x the solution vector.
! w(1) the minimum 2-norm || a.x-b ||.
! istate vector indicating which components of x are active and which are
! at their bounds (see the previous paragraph).
! istate can be supplied to the routine to give it a good starting
! guess for the solution.
! loopA number of iterations taken in the main loop, Loop A.
!-----------------------------------------------------------------------
! Working arrays:
! w dimension n. act dimension m*(mm+2).
! zz dimension m. istate dimension n+1.
!-----------------------------------------------------------------------
! Method: active variable method along the general plan of NNLS by
! Lawson & Hanson, "Solving Least Squares Problems", Prentice-Hall 1974.
! See Algorithm 23.10. Step numbers in comment statements refer to their
! scheme.
! For more details and further uses, see the article
! "Bounded Variable Least Squares: An Algorithm and Applications"
! by Stark and Parker in 1995 Computational Statistics.
!-----------------------------------------------------------------------
! A number of measures are taken to enhance numerical reliability:
! 1. As noted by Lawson and Hanson, roundoff errors in the computation of the
! gradient of the misfit may cause a component on the bounds to appear to
! want to become active, yet when the component is added to the active set,
! it moves away from the feasible region. In this case the component is not
! made active, the gradient of the misfit with respect to a change in that
! component is set to zero, and the program returns to the Kuhn-Tucker test.
! Flag ifrom5 is used in this test, which occurs at the end of Step 6.
! 2. When the least-squares minimizer after Step 6 is infeasible, it is used
! in a convex interpolation with the previous solution to obtain a feasible
! vector. The constant in this interpolation is supposed to put at least
! one component of x on a bound. There can be difficulties:
! 2a. Sometimes, due to roundoff, no interpolated component ends up on
! a bound. The code in Step 11 uses the flag jj, computed in Step 8,
! to ensure that at least the component that determined the
! interpolation constant alpha is moved to the appropriate bound.
! This guarantees that what Lawson and Hanson call `Loop B' is finite.
! 2b. The code in Step 11 also incorporates Lawson and Hanson's feature that
! any components remaining infeasible at this stage (which must be due to
! roundoff) are moved to their nearer bound.
! 3. If the columns of a passed to qr are linearly dependent, the new
! potentially active component is not introduced: the gradient of the
! misfit with respect to that component is set to zero, and control
! returns to the Kuhn-Tucker test.
! 4. When some of the columns of a are approximately linearly dependent,
! we have observed cycling of active components: a component just moved
! to a bound desires immediately to become active again; qr allows it
! to become active and a different component is moved to its bound.
! This component immediately wants to become active, which qr allows, and
! the original component is moved back to its bound. We have taken two
! steps to avoid this problem:
! 4a. First, the column of the matrix a corresponding to the new
! potentially active component is passed to qr as the last column of
! its matrix. This ordering tends to make a component recently moved
! to a bound fail the test mentioned in (1), above.
! 4b. Second, we have incorporated a test that prohibits short cycles.
! If the most recent successful change to the active set was to move
! the component x(jj) to a bound, x(jj) is not permitted to reenter
! the solution at this stage. This test occurs just after checking
! the Kuhn-Tucker conditions, and uses the flag jj, set in Step 8.
! The flag jj is reset after Step 6 if Step 6 was entered from
! Step 5 indicating that a new component has successfully entered the
! active set. The test for resetting jj uses the flag ifrom5,
! which will not equal zero in case Step 6 was entered from Step 5.
! dimension w(n), act(m,min(m,n)+2), zz(m), istate(n+1)
! Local variables
REAL (dp), PARAMETER :: eps = 1.0E-11_dp
INTEGER :: i, iact, ifrom5, it, j, jj, k, k1, kk, ks, mm, mm1, &
nact, nbound, noldb
REAL (dp) :: act(m,m+2), alf, alpha, bad, bdiff, bnorm, bound, &
bsq, obj, resq, ri, sj, w(n), worst, zz(m)
!----------------------First Executable Statement-----------------------
! Step 1. Initialize everything--active and bound sets, initial values, etc.
! Initialize flags, etc.
mm = MIN(m,n)
mm1 = mm + 1
jj = 0
ifrom5 = 0
! Check consistency of given bounds bl, bu.
bdiff = 0.0_dp
DO j = 1, n
bdiff = MAX(bdiff, bu(j)-bl(j))
IF (bl(j) > bu(j)) THEN
WRITE(*, *) ' Inconsistent bounds in BVLS. '
STOP
END IF
END DO
IF (bdiff == 0.0_dp) THEN
WRITE(*, *) ' No free variables in BVLS--check input bounds.'
STOP
END IF
! In a fresh initialization (key = 0) bind all variables at their lower bounds.
! If (key != 0), use the supplied istate vector to initialize the variables.
! istate(n+1) contains the number of bound variables. The absolute values of
! the first nbound=istate(n+1) entries of istate are the indices of the
! bound variables. The sign of each entry determines whether the indicated
! variable is at its upper (positive) or lower (negative) bound.
IF (key == 0) THEN
nbound = n
nact = 0
DO j = 1, nbound
istate(j) = -j
END DO
ELSE
nbound = istate(n+1)
END IF
nact = n - nbound
IF (nact > mm) THEN
WRITE(*, *) ' Too many active variables in BVLS starting solution!'
STOP
END IF
DO k = 1, nbound
j = ABS(istate(k))
IF (istate(k) < 0) x(j) = bl(j)
IF (istate(k) > 0) x(j) = bu(j)
END DO
! In a warm start (key != 0) initialize the active variables to (bl+bu)/2.
! This is needed in case the initial qr results in active variables
! out-of-bounds and Steps 8-11 get executed the first time through.
DO k = nbound + 1, n
kk = istate(k)
x(kk) = (bu(kk)+bl(kk)) / 2
END DO
! Compute bnorm, the norm of the data vector b, for reference.
bsq = SUM( b(1:m)**2 )
bnorm = SQRT(bsq)
!-----------------------------Main Loop---------------------------------
! Initialization complete. Begin major loop (Loop A).
DO loopa = 1, 3 * n
! Step 2.
! Initialize the negative gradient vector w(*).
obj = 0.0_dp
w(1:n) = 0.0_dp
! Compute the residual vector b-a.x , the negative gradient vector
! w(*), and the current objective value obj = || a.x - b ||.
! The residual vector is stored in the mm+1'st column of act(*,*).
DO i = 1, m
ri = b(i) - DOT_PRODUCT( a(i,1:n), x(1:n) )
obj = obj + ri ** 2
w(1:n) = w(1:n) + a(i,1:n) * ri
act(i,mm1) = ri
END DO
! Converged? Stop if the misfit << || b ||, or if all components are
! active (unless this is the first iteration from a warm start).
IF (SQRT(obj) <= bnorm*eps.OR.(loopa > 1.AND.nbound == 0)) THEN
istate(n+1) = nbound
w(1) = SQRT(obj)
RETURN
END IF
! Add the contribution of the active components back into the residual.
DO k = nbound + 1, n
j = istate(k)
act(1:m,mm1) = act(1:m,mm1) + a(1:m,j) * x(j)
END DO
! The first iteration in a warm start requires immediate qr.
IF (loopa == 1 .AND. key /= 0) GO TO 150
! Steps 3, 4.
! Find the bound element that most wants to be active.
120 worst = 0.0_dp
it = 1
DO j = 1, nbound
ks = ABS(istate(j))
bad = w(ks) * SIGN(1,istate(j))
IF (bad < worst) THEN
it = j
worst = bad
iact = ks
END IF
END DO
! Test whether the Kuhn-Tucker condition is met.
IF (worst >= 0.0_dp) THEN
istate(n+1) = nbound
w(1) = SQRT(obj)
RETURN
END IF
! The component x(iact) is the one that most wants to become active.
! If the last successful change in the active set was to move x(iact)
! to a bound, don't let x(iact) in now: set the derivative of the misfit
! with respect to x(iact) to zero and return to the Kuhn-Tucker test.
IF (iact == jj) THEN
w(jj) = 0.0_dp
GO TO 120
END IF
! Step 5.
! Undo the effect of the new (potentially) active variable on the
! residual vector.
IF (istate(it) > 0) bound = bu(iact)
IF (istate(it) < 0) bound = bl(iact)
act(1:m,mm1) = act(1:m,mm1) + bound * a(1:m,iact)
! Set flag ifrom5, indicating that Step 6 was entered from Step 5.
! This forms the basis of a test for instability: the gradient calculation
! shows that x(iact) wants to join the active set; if qr puts x(iact) beyond
! the bound from which it came, the gradient calculation was in error and
! the variable should not have been introduced.
ifrom5 = istate(it)
! Swap the indices (in istate) of the new active variable and the
! rightmost bound variable; `unbind' that location by decrementing nbound.
istate(it) = istate(nbound)
nbound = nbound - 1
nact = nact + 1
istate(nbound+1) = iact
IF (mm < nact) THEN
WRITE(*, *) ' Too many free variables in BVLS!'
STOP
END IF
! Step 6.
! Load array act with the appropriate columns of a for qr. For added
! stability, reverse the column ordering so that the most recent addition to
! the active set is in the last column. Also copy the residual vector from
! act(., mm1) into act(., mm1+1).
150 DO i = 1, m
act(i,mm1+1) = act(i,mm1)
DO k = nbound + 1, n
j = istate(k)
act(i, nact+1-k+nbound) = a(i,j)
END DO
END DO
CALL qr(m, nact, act, act(:,mm1+1), zz, resq)
! Test for linear dependence in qr, and for an instability that moves the
! variable just introduced away from the feasible region (rather than into
! the region or all the way through it).
! In either case, remove the latest vector introduced from the active set
! and adjust the residual vector accordingly.
! Set the gradient component (w(iact)) to zero and return to the Kuhn-Tucker
! test.
IF (resq < 0.0_dp .OR. (ifrom5 > 0 .AND. zz(nact) > bu(iact)) .OR. &
(ifrom5 < 0 .AND. zz(nact) < bl(iact))) THEN
nbound = nbound + 1
istate(nbound) = istate(nbound) * SIGN(1.0D0, x(iact)-bu(iact))
nact = nact - 1
act(1:m,mm1) = act(1:m,mm1) - x(iact) * a(1:m,iact)
ifrom5 = 0
w(iact) = 0.0_dp
GO TO 120
END IF
! If Step 6 was entered from Step 5 and we are here, a new variable
! has been successfully introduced into the active set; the last
! variable that was fixed at a bound is again permitted to become active.
IF (ifrom5 /= 0) jj = 0
ifrom5 = 0
! Step 7. Check for strict feasibility of the new qr solution.
DO k = 1, nact
k1 = k
j = istate(k+nbound)
IF (zz(nact+1-k) < bl(j) .OR. zz(nact+1-k) > bu(j)) GO TO 210
END DO
DO k = 1, nact
j = istate(k+nbound)
x(j) = zz(nact+1-k)
END DO
! New iterate is feasible; back to the top.
CYCLE
! Steps 8, 9.
210 alpha = 2.0_dp
alf = alpha
DO k = k1, nact
j = istate(k+nbound)
IF (zz(nact+1-k) > bu(j)) alf = (bu(j)-x(j)) / ( zz(nact+1-k)-x(j))
IF (zz(nact+1-k) < bl(j)) alf = (bl(j)-x(j)) / ( zz(nact+1-k)-x(j))
IF (alf < alpha) THEN
alpha = alf
jj = j
sj = SIGN(1.0_dp, zz(nact+1-k)-bl(j))
END IF
END DO
! Step 10
DO k = 1, nact
j = istate(k+nbound)
x(j) = x(j) + alpha * (zz(nact+1-k)-x(j))
END DO
! Step 11.
! Move the variable that determined alpha to the appropriate bound.
! (jj is its index; sj is + if zz(jj)> bu(jj), - if zz(jj) 0.0_dp)) THEN
! Move x(j) to its upper bound.
x(j) = bu(j)
istate(k+noldb) = istate(nbound+1)
istate(nbound+1) = j
nbound = nbound + 1
act(1:m,mm1) = act(1:m,mm1) - bu(j) * a(1:m,j)
ELSE IF (x(j)-bl(j) <= 0.0_dp .OR. (j == jj .AND. sj < 0.0_dp)) THEN
! Move x(j) to its lower bound.
x(j) = bl(j)
istate(k+noldb) = istate(nbound+1)
istate(nbound+1) = -j
nbound = nbound + 1
act(1:m,mm1) = act(1:m,mm1) - bl(j) * a(1:m,j)
END IF
END DO
nact = n - nbound
! If there are still active variables left repeat the qr; if not,
! go back to the top.
IF (nact > 0) GO TO 150
END DO
WRITE(*, *) ' BVLS fails to converge! '
STOP
END SUBROUTINE bvls
!======================================================================
SUBROUTINE qr(m, n, a, b, x, resq)
!======================================================================
INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: a(m,n)
REAL (dp), INTENT(IN OUT) :: b(m)
REAL (dp), INTENT(OUT) :: x(n)
REAL (dp), INTENT(OUT) :: resq
!$$$$ calls no other routines
! Relies on FORTRAN77 do-loop conventions!
! Solves over-determined least-squares problem ax ~ b
! where a is an m by n matrix, b is an m-vector.
! resq is the sum of squared residuals of optimal solution. Also used
! to signal error conditions - if -2 , system is underdetermined, if
! -1, system is singular.
! Method - successive Householder rotations. See Lawson & Hanson -
! Solving Least Squares Problems (1974).
! Routine will also work when m=n.
!***** CAUTION - a and b are overwritten by this routine.
! Local variables
REAL (dp) :: const, total, dot, qv1, sq, u1
INTEGER :: i, j, j1, jj
resq = -2.0
IF (m < n) RETURN
resq = -1.0
! Loop ending on 1800 rotates a into upper triangular form.
DO j = 1, n
! Find constants for rotation and diagonal entry.
sq = SUM( a(j:m,j) ** 2 )
IF (sq == 0.0_dp) RETURN
qv1 = -SIGN(SQRT(sq), a(j,j))
u1 = a(j,j) - qv1
a(j,j) = qv1
j1 = j + 1
! Rotate remaining columns of sub-matrix.
DO jj = j1, n
dot = u1 * a(j,jj) + DOT_PRODUCT( a(j1:m,jj), a(j1:m,j) )
const = dot / ABS(qv1*u1)
a(j1:m,jj) = a(j1:m,jj) - const * a(j1:m,j)
a(j,jj) = a(j,jj) - const * u1
END DO
! Rotate b vector.
dot = u1 * b(j) + DOT_PRODUCT( a(j1:m,j), b(j1:m) )
const = dot / ABS(qv1*u1)
b(j) = b(j) - const * u1
b(j1:m) = b(j1:m) - const * a(j1:m,j)
END DO
! Solve triangular system by back-substitution.
DO i = n, 1, -1
total = b(i) - DOT_PRODUCT( a(i,i+1:n), x(i+1:n) )
IF (a(i,i) == 0.0_dp) RETURN
x(i) = total / a(i,i)
END DO
! Find residual in overdetermined case.
resq = SUM( b(n+1:m) ** 2 )
RETURN
END SUBROUTINE qr
!______________________________________________________________________
END MODULE BoundedLeastSquares
PROGRAM Test_BVLS
USE BoundedLeastSquares
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
INTEGER, PARAMETER :: ncases = 100, ncols = 10
REAL (dp) :: a(ncases,ncols), y(ncases), bl(ncols), bu(ncols), &
beta(ncols), x(ncols), e
INTEGER :: case, istate(ncols+1), j, key, loopa
! Generate artificial data satisfying
! Y = A.beta + noise
CALL RANDOM_NUMBER(beta)
beta = 4.0*(beta - 0.5)
CALL RANDOM_NUMBER(a)
DO case = 1, ncases
CALL RANDOM_NUMBER(e)
y(case) = DOT_PRODUCT( a(case, :), beta ) + 0.1*(e - 0.5)
END DO
key = 0
bl = 0.0_dp
bu = 1.0_dp
CALL bvls(key, ncases, ncols, a, y, bl, bu, x, istate, loopa)
WRITE(*, *) ' Column Original beta Solution'
DO j = 1, ncols
WRITE(*, '(i5, 2f14.3)') j, beta(j), x(j)
END DO
WRITE(*, '(a, i4)') ' No. of iterations = ', loopa
WRITE(*, *)
STOP
END PROGRAM Test_BVLS