! Compliments of netlib Sun Jul 6 04:10:36 EDT 1986 ! algorithm 615, collected algorithms from acm. ! algorithm appeared in acm-trans. math. software, vol.10, no. 2, june, 1984, ! pp. 202-206. Authors: Armstrong, R.D., Beck, P.O. & Kung, M.T. MODULE toms615 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) ! Replace the COMMON statements in TOMS615 REAL (dp), SAVE :: acu, big, bsave(210), lu(20,20), pie(210), sad, sadt, & savtot(210), sig(6000), sigma(300), zlow, zsave(20) INTEGER, SAVE :: ibase(20), ik, ik1, incol(20), index(20), inprob(210), & ipar(20), label(210), level LOGICAL, SAVE :: direct, intl CONTAINS SUBROUTINE kbest (x, y, m, n, iter, ifault, popt, minin, nmax, mmax, bval, & idex, istat, zl) !*********************************************************************** ! The purpose of this program is to determine the best subset of parameters ! to fit a linear regression under an least absolute value criterion. ! This program utilizes the simplex method of linear programming within a ! branch-and-bound algorithm to solve the best subset problem. ! The algorithm is based on the publication: ! Armstrong,R.D. and M.T. Kung "An algorithm to select the best subset for ! a least absolute value regression problem", Optimization in Statistics, ! TIMS Studies of the Management Sciences, vol.33, 931-936 (1982). ! formal parameters ! x real array input: values of independent variables such that ! (nmax,mmax) each row corresponds to an observation ! y real array input: values of the dependent variables ! (nmax) ! m integer input: number of dependent variables ! n integer input: number of observations ! iter integer output: number of iterations ! ifault integer output: failure indicator ! = 0 normal termination ! = 1 observation matrix does not have full ! row rank (rank m) ! = 2 problem size out of range ! = 3 no pivot element found implying near ! singular basis ! popt real input: percentage deviation from ! optimality allowed ! minin integer input: minimum number of parameters in the ! model. best subset of size minin t ! m is obtained. ! nmax integer input: dimension of rows in x (also y) ! mmax integer input: dimension of columns in x ! bval real array output: array of optimal beta values for each ! subset. The beta values for the subset ! of size m are stored in positions ! bval(1), bval(2), ..., bval(m), for the ! subset of size m-1 the values are stored ! in positions bval(m+1), bval(m+2), ..., ! bval(2m-1). In general, the beta values ! for the optimal subset of size k are ! stored in positions bval(l), ..., ! bval(l-k+1) where ! l = (m*(m+1)-k*(k+1))/2 + 1 ! idex integer array output: beta index set for the optimal subset. ! (((mmax+1)*mmax)/2) This array is a parallel array for bval; ! i.e., if bval(j) = 2. and idex(j) = 7 ! then beta(7) = 2.7 in the associated ! optimal subset. ! istat integer array input: parameter status array. ! (mmax) ! 1 if beta(j) is required ! in every model ! istat(j) = ! 0 otherwise ! zl real array output: best objective value for each subse ! (mmax) ! zl(j) gives the best objective valu ! for the subset with m-j+1 parameter !*********************************************************************** ! implementation notes: ! 1. the routine uses two machine dependent values acu and big. ! (these are set at the beginning of this subroutine) ! a) acu is used to test for "zero". acu should be set to approximately ! 100 * the relative machine accuracy of the system in use. ! b) big is used to initialize the array zl (described above). ! big should be set to the largest floating point value assignable. ! 2. both single and double precision versions are supplied. This version ! is in double precision. ! 3. array dimensions ! the code is currently dimensioned to solve problems with up to 20 ! parameters and 300 observations. the dimension sizes are ! determined as follows ! 20 mmax ! 300 nmax ! 210 (nmax+1)*mmax/2 ! 6000 nmax*mmax ! 10 maximum value of ik . ! if the dimension sizes are changed the common and dimension ! statements will need to be modified in each subroutine. !******************************************************************** ! subroutines: ! calbet - back-solves a sysytem of equations ! calcpi - forward-solves a sysytem of equations ! kbest - the driver ! l1norm - solves the initial regression problem with all ! parameters included in the model ! phase2 - solves the current regression problem with a primal algorithm ! setup - determines the form of the current problem by ! choosing the parameter to leave based on a penalty ! update - updates lu decomposition matrix !*********************************************************************** ! description of variables: ! ik: length of candidate list ! beta: estimates to the current subproblem ! sad: the minimum total absolute deviation ! ibase: the index array of columns of the basis ! n,m,nmax,mmax,x,y, are not changed in the subroutine; ! sad is updated at each iteration ! lu: the lu decomposition of the current basis ! index: the index array of rows of lu ! tot: the current rhs of the dual problem ! sigma: indicator array to specify whether a nonbasic dual variable ! is at upper or lower bound(+1 implies upper; -1 implies lower) ! inext: a local array for the sort routine ! rhs: a local array for the calbet routine ! ipar(i) = -k if the k-th parameter (beta) is out of model at level i ! = k if the k-th parameter is in ! label(i) saves the indices of the basic observations in the ! predecessor path ! zsave(i) saves the objective values on the predecessor path !********************************************************************** ! Conversion to be ELF90 compatible by Alan Miller ! e-mail: amiller @ bigpond.net.au ! http://users.bigpond.net.au ! Latest revision - 8 December 1998 IMPLICIT NONE REAL (dp), INTENT(IN) :: x(300,20), y(300), popt INTEGER, INTENT(IN) :: m, n, nmax, mmax INTEGER, INTENT(IN OUT) :: minin, istat(20) INTEGER, INTENT(OUT) :: iter, ifault, idex(210) REAL (dp), INTENT(OUT) :: bval(210), zl(20) ! Local variables REAL (dp) :: beta(20), tot(20), pi(20), repp(20) INTEGER :: i, j, jinm, k2, k3, kkk1, l, nfree, numin REAL (dp) :: popt1, zero = 0.0_dp ! assign values to machine dependent constants ! acu = SQRT(d1mach(4)) ! big = d1mach(2) acu = SQRT( 2.0_dp * EPSILON(1.0_dp) ) big = HUGE(1.0_dp) ifault = 0 ik = 5 ik1 = ik-1 popt1 = (100. - popt)/100. ! test problem size IF (n <= nmax .AND. m <= mmax) GO TO 10 ifault = 2 RETURN 10 DO i = 1,m index(i) = i tot(i) = zero incol(i) = i END DO CALL l1norm (x, y, n, m, ifault, iter, beta, pi, tot) IF (ifault > 0) RETURN ! initialization ! save the initial solution jinm = 0 DO i = 1,m bsave(i) = beta(i) bval(i) = beta(i) pie(i) = pi(i) label(i) = ibase(i) idex(i) = i inprob(i) = i savtot(i) = tot(i) zl(i) = big jinm = jinm+istat(i) END DO sig(1:n) = sigma(1:n) level = 0 direct = .true. IF (jinm > minin) minin = jinm ! initialize number of free parameters in model nfree = m-jinm zl(m) = sad zsave(m) = sad IF (minin == m) RETURN ! numin gives the number of parameters in the model numin = m ! start the main loop 50 numin = numin-1 level = level+1 IF (numin < minin .OR. nfree == 0) GO TO 120 k2 = (m*(m+1) - (numin+1)*(numin+2))/2 kkk1 = numin+1 zlow = zl(numin)*popt1 IF (direct) GO TO 80 DO i = 1,kkk1 j = k2+i index(i) = i ibase(i) = label(j) incol(i) = inprob(j) pi(i) = pie(j) beta(i) = bsave(j) tot(i) = savtot(j) END DO ! load in the bounds from the immediate predecessor k3 = (m-kkk1)*n DO i = 1,n k3 = k3+1 sigma(i) = sig(k3) END DO sad = zsave(kkk1) ! set up a new problem 80 CALL setup (x, kkk1, ifault, istat, beta, tot, pi, repp) IF (ifault > 0) RETURN nfree = nfree - 1 ! call primal l.p. code CALL phase2 (x, y, numin, n, iter, ifault, beta, tot, pi, repp) ! save the solution data for later recall k2 = k2 + kkk1 + 1 j = k2 DO i = 1,numin l = k2+i pie(l) = pi(i) label(l) = ibase(i) bsave(l) = beta(i) inprob(l) = incol(i) savtot(l) = tot(i) END DO ! save the objective value zsave(numin) = sad ! save the nonbasic bounds k3 = (m-numin)*n DO i = 1,n k3 = k3+1 sig(k3) = sigma(i) END DO direct = .true. ! check for a better solution IF (sadt >= zlow) GO TO 50 ! save the better solution zl(numin) = sad k2 = j DO i = 1,numin l = k2 + i idex(l) = incol(i) bval(l) = beta(i) END DO ! go to top of loop GO TO 50 120 numin = numin + 1 direct = .false. ! must work back up the tree level = level - 1 130 IF (ipar(level) > 0) GO TO 140 ipar(level) = -ipar(level) j = ipar(level) istat(j) = 1 numin = numin + 1 GO TO 50 140 j = ipar(level) istat(j) = 0 nfree = nfree + 1 level = level - 1 IF (level > 0) GO TO 130 RETURN END SUBROUTINE kbest SUBROUTINE setup (x, m, ifault, istat, beta, tot, pi, repp) ! this subroutine determines the form of the current subproblem IMPLICIT NONE REAL (dp), INTENT(IN) :: x(300,20), beta(20) REAL (dp), INTENT(IN OUT) :: tot(20), pi(20), repp(20) INTEGER, INTENT(OUT) :: ifault INTEGER, INTENT(IN OUT) :: m, istat(20) ! Local variables INTEGER :: i, ifixi, ii, iout, isave, kkk, l, leave, ll REAL (dp) :: beside, one = 1.0_dp, penalt, ratio, rho, rhs(20), side, srepp, & test, zero = 0.0_dp ! if the immediate predecessor is in memory go on IF (direct) GO TO 10 ! reconstruct the lu decomposition kkk = 1 CALL update (kkk, ifault, x, 1, m) IF (ifault > 0) RETURN ! recalculate the basic pi values kkk = 0 CALL calcpi (kkk, pi, tot, m) ! determine parameter to leave based on penalty 10 penalt = big kkk = 0 rhs(1:m) = zero DO i = 1,m ii = index(i) ll = incol(ii) IF (istat(ll) == 1) CYCLE rho = SIGN(1._dp, beta(ii)) rhs(ii) = one kkk = i-1 CALL calcpi (kkk, repp, rhs, m) rhs(ii) = zero test = big DO l = 1,m IF (ABS(repp(l)) <= acu) CYCLE srepp = SIGN(1._dp,repp(l)) ratio = ABS((1._dp - rho*srepp*pi(l))/repp(l)) IF (ratio >= test) CYCLE side = rho*srepp test = ratio iout = l END DO ! see if the penalty is less IF (test*ABS(beta(ii)) >= penalt) CYCLE leave = iout ifixi = ii beside = side penalt = test*ABS(beta(ii)) END DO ! update the objective value sad = sad+penalt ! parameter incol(ifixi) will leave the model ! pi(ibase(leave)) will leave the basis and go to bound ! switch unwanted parameter and pi to the end of list isave = incol(ifixi) incol(ifixi) = incol(m) incol(m) = isave tot(ifixi) = tot(m) ! label this node and flag the outgoing parameter ipar(level) = -isave istat(isave) = -1 isave = ibase(leave) ibase(leave) = ibase(m) ibase(m) = isave ! place the outgoing pi at the correct bound sigma(isave) = beside ! form a new basis m = m-1 DO i = 1,m index(i) = i END DO kkk = 1 CALL update (kkk, ifault, x, 1, m) IF (ifault > 0) RETURN ! update the tot array DO i = 1,m ll = incol(i) tot(i) = tot(i) - beside*x(isave,ll) END DO RETURN END SUBROUTINE setup SUBROUTINE phase2 (x, y, m, n, iter, ifault, beta, tot, pi, repp) ! solves lp with primal algorithm IMPLICIT NONE INTEGER, INTENT(IN) :: m, n INTEGER, INTENT(IN OUT) :: iter INTEGER, INTENT(OUT) :: ifault REAL (dp), INTENT(IN) :: x(300,20), y(300) REAL (dp), INTENT(IN OUT) :: tot(20) REAL (dp), INTENT(OUT) :: beta(20), pi(20), repp(20) ! Local variables INTEGER :: i, iin, iout, j, k, kan(10), kkk, krow, l, ll REAL (dp) :: one = 1.0_dp, ratio, rho, rhs(20), srepp, test, xtwo, & z(10), zero = 0.0_dp, zc, zz ! z(i) stores the reduced cost values on the list ! kan(i) stores the row number on the list ! ik is the length of the list ! calculate the simplex multipliers DO i = 1,m k = ibase(i) pi(i) = y(k) END DO kkk = 0 CALL calbet (kkk, beta, pi, m) ! recalculate the basic pi values kkk = 0 CALL calcpi (kkk, pi, tot, m) ! store the objective value used for termination check sadt = sad IF (sad >= zlow) RETURN ! generate the candidate list 20 z(1:ik) = acu ! calculate the reduced costs DO j = 1,n IF (sigma(j) == zero) CYCLE zc = zero DO i = 1,m ll = incol(i) zc = zc+beta(i)*x(j,ll) END DO zz = (zc - y(j))*sigma(j) IF (zz <= z(ik)) CYCLE ! rank the values in z(i), in descending order DO k = 1,ik1 IF (zz <= z(k)) CYCLE i = ik 50 z(i) = z(i-1) kan(i) = kan(i-1) i = i-1 IF (i > k) GO TO 50 z(k) = zz kan(k) = j CYCLE END DO z(ik) = zz kan(ik) = j END DO IF (z(1) <= acu) RETURN krow = 1 zz = z(1) iin = kan(1) zc = zz*sigma(iin) GO TO 100 80 krow = krow+1 IF (krow > ik) GO TO 20 IF (z(krow) <= acu) GO TO 20 ! determine the entering variable from the candidate list iin = kan(krow) zc = zero DO i = 1,m ll = incol(i) zc = zc + beta(i)*x(iin,ll) END DO zc = zc - y(iin) zz = zc*sigma(iin) IF (zz <= acu) GO TO 80 100 rho = sigma(iin) ! rho: sign of the incoming reduced cost ! iin: incoming candidate ! find the representation of the entering variable DO i = 1,m ll = incol(i) rhs(i) = x(iin,ll) END DO kkk = 0 CALL calcpi (kkk, repp, rhs, m) ! calculate the min ratio test to determine the leaving variable test = 2.0 DO i = 1,m IF (ABS(repp(i)) <= acu) CYCLE srepp = SIGN(1.0_dp,repp(i)) ratio = ABS((1._dp - rho*srepp*pi(i))/repp(i)) IF (ratio >= test) CYCLE test = ratio iout = i END DO ! perform fathoming test before pivot sadt = sad + test*zz IF (sadt >= zlow) RETURN sad = sadt pi(1:m) = pi(1:m) + test*repp(1:m)*rho IF (test < 2.0) GO TO 150 sigma(iin) = -sigma(iin) xtwo = sigma(iin) + sigma(iin) DO i = 1,m ll = incol(i) tot(i) = tot(i) - xtwo*x(iin,ll) END DO GO TO 80 150 k = ibase(iout) sigma(k) = SIGN(1.0_dp, pi(iout)) pi(iout) = sigma(iin)*(1.0-test) ibase(iout) = iin DO i = 1,m ll = incol(i) tot(i) = tot(i) + sigma(iin)*x(iin,ll) - sigma(k)*x(k,ll) END DO kkk = iout CALL update (iout, ifault, x, 1, m) rhs(1:m) = zero rhs(kkk) = one kkk = kkk - 1 CALL calbet (kkk, repp, rhs, m) DO l = 1,m beta(l) = beta(l) - zc*repp(l) END DO sigma(iin) = zero iter = iter + 1 GO TO 80 END SUBROUTINE phase2 SUBROUTINE l1norm (x, y, n, m, ifault, iter, beta, pi, tot) ! this routine solves the initial regression problem with all ! the parameters included in the model IMPLICIT NONE INTEGER, INTENT(IN) :: n, m INTEGER, INTENT(OUT) :: ifault, iter REAL (dp), INTENT(IN) :: x(300,20), y(300) REAL (dp), INTENT(IN OUT) :: tot(20) REAL (dp), INTENT(OUT) :: beta(20), pi(20) ! Local variables REAL (dp) :: d(300), delta(300), price(300), rhs(20) INTEGER :: i, iin, inext(300), iout, ipt, j, k, k1, kkk, kount, l REAL (dp) :: ahalf, aone, one = 1.0_dp, ratio, rho, subt, t, test, total, & val, zero = 0.0_dp, zzz ! initial settings iter = 0 aone = one + acu ahalf = .5 + acu ! set up initial lu decomposition intl = .true. k = 1 CALL update (k, ifault, x, n, m) IF (ifault /= 0) RETURN intl = .false. DO i = 1,m k1 = ibase(i) rhs(i) = y(k1) END DO kkk = 0 CALL calbet (kkk, beta, rhs, m) ! calculate initial d, tot and sigma vectors DO j = 1,n val = DOT_PRODUCT( beta(1:m), x(j,1:m) ) d(j) = y(j) - val sigma(j) = SIGN(1._dp, d(j)) END DO DO j = 1,m rhs(j) = zero kkk = ibase(j) sigma(kkk) = zero END DO DO j = 1,n DO i = 1,m tot(i) = tot(i) - sigma(j)*x(j,i) END DO END DO ! main iterative loop begins 70 kkk = 0 CALL calcpi (kkk, pi, tot, m) t = aone k = 0 DO j = 1,m IF (ABS(pi(j)) < t) CYCLE k = j t = ABS(pi(j)) rho = -SIGN(1._dp, pi(j)) END DO IF (k == 0) GO TO 220 kkk = k - 1 rhs(k) = one CALL calbet (kkk, beta, rhs, m) rhs(k) = zero DO i = 1,n delta(i) = zero IF (sigma(i) == zero) CYCLE DO j = 1,m delta(i) = delta(i) + beta(j)*x(i,j) END DO delta(i) = rho*delta(i) END DO ! perform partial sort of ratios t = t*.5 kount = 0 ratio = big total = ahalf subt = zero DO i = 1,n IF (delta(i)*sigma(i) <= acu) CYCLE test = d(i)/delta(i) IF (test >= ratio) CYCLE total = total + ABS(delta(i)) IF (total - subt >= t) GO TO 110 ! insert i in list kount = kount+1 price(kount) = test inext(kount) = i CYCLE ! update sum and kick iin out of the list 110 total = total - subt ratio = test ipt = 0 kkk = 0 ! identify a new iin 120 kkk = kkk + 1 IF (kkk > kount) GO TO 130 IF (price(kkk) <= ratio) GO TO 120 ratio = price(kkk) ipt = kkk GO TO 120 130 IF (ipt == 0) GO TO 150 ! switch values kkk = inext(ipt) subt = ABS(delta(kkk)) IF (total - subt < t) GO TO 140 price(ipt) = price(kount) inext(ipt) = inext(kount) kount = kount - 1 GO TO 110 140 iin = inext(ipt) inext(ipt) = i price(ipt) = test CYCLE 150 iin = i subt = ABS(delta(i)) END DO ! update basic indicators DO j = 1,kount kkk = inext(j) zzz = sigma(kkk) DO l = 1,m subt = zzz*x(kkk,l) tot(l) = tot(l) + subt + subt END DO sigma(kkk) = -sigma(kkk) END DO iout = ibase(k) delta(iout) = rho DO l = 1,m tot(l) = tot(l) + rho*x(iout,l) + sigma(iin)*x(iin,l) END DO sigma(iout) = -rho ibase(k) = iin CALL update (k, ifault, x, 1, m) d(1:n) = d(1:n) - ratio*delta(1:n) sigma(iin) = zero iter = iter + 1 GO TO 70 ! calculate optimal beta and sum of absolute deviations 220 DO i = 1,m k1 = ibase(i) rhs(i) = y(k1) END DO kkk = 0 CALL calbet (kkk, beta, rhs, m) sad = SUM( ABS(d(1:n)) ) RETURN END SUBROUTINE l1norm SUBROUTINE calbet (kkk, beta, rhs, m) ! subroutine calbet for back-solving system of equations IMPLICIT NONE INTEGER, INTENT(IN) :: m INTEGER, INTENT(IN OUT) :: kkk REAL (dp), INTENT(IN) :: rhs(20) REAL (dp), INTENT(OUT) :: beta(20) ! Local variables INTEGER :: i, ii, ii1, k, k1, k2, kk, kkk1, m1 REAL (dp) :: zero = 0.0_dp IF (m > 1) GO TO 10 k = index(1) beta(k) = rhs(1)/lu(k,1) RETURN 10 m1 = m-1 IF (kkk == 0) GO TO 30 DO i = 1,kkk k = index(i) beta(k) = zero END DO 30 kkk = kkk + 1 k = index(kkk) beta(k) = rhs(kkk)/lu(k,kkk) IF (kkk == m) GO TO 60 kkk1 = kkk + 1 DO ii = kkk1,m k = index(ii) beta(k) = rhs(ii) ii1 = ii - 1 DO i = kkk,ii1 kk = index(i) beta(k) = beta(k) - lu(kk,ii)*beta(kk) END DO beta(k) = beta(k)/lu(k,ii) END DO 60 DO ii = 1,m1 k1 = m - ii k = index(k1) DO i = 1,ii kk = m - i + 1 k2 = index(kk) beta(k) = beta(k) - lu(k2,k1)*beta(k2) END DO END DO RETURN END SUBROUTINE calbet SUBROUTINE update (kkk, ifault, x, n, m) ! subroutine update for updating lu decomposition matrix IMPLICIT NONE INTEGER, INTENT(IN) :: n, m INTEGER, INTENT(IN OUT) :: kkk INTEGER, INTENT(OUT) :: ifault REAL (dp), INTENT(IN) :: x(300,20) INTEGER :: i, icol, ii, ii1, irow, isave, j, k, kk, ll REAL (dp) :: pivot, subt irow = 0 DO ii = kkk,m IF (intl) GO TO 10 irow = ibase(ii) GO TO 20 10 irow = irow + 1 ibase(ii) = irow IF (irow > n) ifault = 1 IF (ifault /= 0) RETURN 20 DO i = 1,m ll = incol(i) lu(i,ii) = x(irow,ll) END DO ! set up representation of incoming row IF (ii == 1) GO TO 60 ii1 = ii - 1 DO icol = 1,ii1 k = index(icol) subt = lu(k,ii) j = icol + 1 DO i = j,m k = index(i) lu(k,ii) = lu(k,ii) - subt*lu(k,icol) END DO END DO ! find maximum entry 60 pivot = acu kk = 0 DO i = ii,m k = index(i) IF (ABS(lu(k,ii)) <= pivot) CYCLE pivot = ABS(lu(k,ii)) kk = i END DO IF (kk == 0) GO TO 10 ! switch order isave = index(kk) index(kk) = index(ii) index(ii) = isave ! put in columns of lu one at a time IF (ii == m) CYCLE j = ii + 1 DO i = j,m k = index(i) lu(k,ii) = lu(k,ii)/lu(isave,ii) END DO END DO kkk = irow RETURN END SUBROUTINE update SUBROUTINE calcpi (kkk, pi, tot, m) ! this routine forward solves a linear system IMPLICIT NONE INTEGER, INTENT(IN) :: m INTEGER, INTENT(IN OUT) :: kkk REAL (dp), INTENT(IN) :: tot(20) REAL (dp), INTENT(OUT) :: pi(20) ! Local variables INTEGER :: i, ii, ii1, k, k1, k2, kkk1, m1 REAL (dp) :: zero = 0.0_dp m1 = m-1 IF (m > 1) GO TO 10 k = index(m) pi(m) = tot(k)/lu(k,m) RETURN 10 pi(1:kkk) = zero kkk = kkk + 1 k = index(kkk) pi(kkk) = tot(k) IF (kkk == m) GO TO 60 kkk1 = kkk + 1 DO ii = kkk1, m k = index(ii) pi(ii) = tot(k) ii1 = ii - 1 pi(ii) = pi(ii) - DOT_PRODUCT( lu(k, kkk:ii1), pi(kkk:ii1) ) END DO 60 k = index(m) pi(m) = pi(m)/lu(k,m) DO ii = 1,m1 k1 = m - ii k = index(k1) DO i = 1,ii k2 = m - i + 1 pi(k1) = pi(k1) - lu(k,k2)*pi(k2) END DO pi(k1) = pi(k1)/lu(k,k1) END DO RETURN END SUBROUTINE calcpi END MODULE toms615