 # document.write(pagetitle) Home Download Tutorial 1 Tutorial 2 Tutorial 3 Tutorial 4 Tutorial 5
```;+
; NAME:              BVLS_PS
;
; PURPOSE:           Bounded value least-square minimization (BVLS)
;
; USAGE:
;                    bvls_ps, a, b, bnd, x, istate, loopa,
;                       RNORM=rnorm
;                       ITMAX=itmax,
;                       EPS=eps
;                       STATUS=status, /QUIET, /DOUBLE,
;
; ARGUMENTS:
;       A   [input]
;          On entry A contains the M by N matrix, A.
;       B   [input]
;          On entry B contains the M-vector, B.
;       BND [input]
;          BND[0,J] is the lower bound for X[J].
;          BND[1,J] is the upper bound for X[J].
;          Require:  BND[0,J]  <=  BND[1,J].
;       X   [output]
;          Solution N-vector.
;       ISTATE [optional input/output]
;          (N+1)-vector containing the state of each variable. If it is
;          given in input, it indicates a "warm" start, i.e. the program
;          will start from a guess of which variables are active (i.e.
;          strictly within their bounds). In output, it contains the final
;          state.
;          This may be used when a series of related problems is treated
;          (pass the vector computed at a preious call for the next).
;          istate[n] is the total number of components at their bounds.
;          The absolute values of the first nbound entries of  istate  are the
;          indices of these `bound' components of  x. The sign of istate[j]
;          j=0...nbound-1, indicates whether  x[|istate[j]|-1] 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,...,n-1  contain the indices of the active
;          components of x.
;
; KEYWORDS:   [Optional input]
;       LOOPA=loopa
;          number of iterations taken in the main loop, Loop A.
;       ITMAX=itmax,
;          Maximum number of iterations permitted. Defaulted to 3*N, this is
;          usually larger than required.
;       EPS=eps,
;          Determines the relative linear dependence of a column vector for a
;          variable moved from its initial value.  The default value is
;          EPS=(machar()).eps (single precision).
;       /DOUBLE,
;          Setting this keyword forces some computations to be performed in
;          double precision. This may be required if: (i) N is large (>10^6)
;          or (ii) if the A or B values span a large amplitude (e.g. 10^20).
;          This slows the computation, and is not required in most cases.
;       /QUIET
;          Verbosity control. By default the procedure prints the messages.
;
; KEYWORDS:   [Optional output]
;       RNORM=rnorm,
;         Euclidean norm of the residual vector, || b - A*X ||.
;       STATUS=status,
;         Indicates status on return.
;            0   Solution completed.
;            1   A is not a 2D array
;            2   B is not a vector
;            3   Sizes of A and B and not compatible
;            4   BND is not a 2D array
;            5   Size of BND is not consistent with B
;            6   Some lower bounds are excessively large
;            7   Some upper bounds are excessively small
;            8   Inconsistent bounds
;            9   No free variable
;           10   Too many active variables in input STATE
;           11   Too many free variables
;           12   Failed to converge (reached ITMAX)
;
; DESCRIPTION:
;   BVLS_PS solves linear least-squares problems with upper and lower bounds
;   on the variables, using an active set strategy. It 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.
;
;   The unconstrained least-squares problems for each candidate set of free
;   variables are solved using the QR decomposition. BVLS has a "warm-start"
;   feature permitting some of the variables to be initialized at their upper
;   or lower bounds, which speeds the solution of a sequence of related
;   problems.
;
;   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.
;
; AUTHOR:
;   Robert L. Parker and Philip B. Stark   (Version 3/19/90)
;   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.
;
; HISTORY:
;   The original Fortran77 code, can be found in the authors web
;   pages: http://www.stat.berkeley.edu/~stark/
;
;   It was converted to Fortran90 by Alan Miller (2000-11-23  Time: 23:41:42)
;   This can be found from: http://jblevins.org/mirror/amiller/
;
;   The GDL/IDL transcoding of this latter was completed by Philippe Prugniel
;   on 2011 January 18th. It use LA_LEAST_SQUARES (IDL) to perform the QR
;   factorization.
;   This routine is used in ULySS for replace the Lawson and Hanson version,
;   BVLS_LH. It is about 3 times faster, and using the warm start capability
;   provides another factor 3 gain when a minimization with several components
;-
; CATEGORY:    ULY_UTIL
;------------------------------------------------------------------------------

;--------------------Bounded Variable Least Squares---------------------
;
;  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
;
;  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.

pro bvls_ps, a, b, bnd, x, istate, \$
RNORM=rnorm,    \$
DOUBLE=double,  \$
EPS=eps,        \$
ITMAX=itmax,    \$
LOOP=loopa,     \$
STATUS=status,  \$
QUIET=quiet

;----------------------First Executable Statement-----------------------
;  Step 1.  Check the dimension of the input arrays, initial values, etc.
;  Initialize everything--active and bound sets

status = 0
if n_elements(eps) eq 0 then eps = (machar()).eps

;  Even for the double precision computations, avoid some unreasonable values
huge = (machar()).xmax             ;  single precision largest usable number

sz = size(a)
if sz ne 2 then begin
status = 1
if not keyword_set(quiet) then message, /INFO, 'Matrix A should be a 2D array'
return
endif
m = sz
n = sz
sz = size(b)
if sz ne 1 then begin
status = 2
if not keyword_set(quiet) then message, /INFO, 'B should be a vector'
return
endif
if sz ne m then begin
status = 3
if not keyword_set(quiet) then \$
message, /INFO, 'The length of B ('+strtrim(sz,2)+ \$
') is not consitent with A ('+strtrim(m,2)+')'
return
endif
sz = size(bnd)
if sz ne 2 then begin
status = 4
if not keyword_set(quiet) then message, /INFO, 'BND should be 2D array'
return
endif
if sz ne 2 then begin
status = 5
if not keyword_set(quiet) then message, /INFO, 'BND should be a 2xN array'
return
endif
if sz ne n then  begin
status = 5
if not keyword_set(quiet) then message, /INFO, 'BND should be a 2xN array'
return
endif

if n_elements(itmax) eq 0 then itmax = 3 * n

bl = bnd[0,*]
nl = where(finite(bl) eq 0, cnt)  ; Take non-finite values as 'No bound'
if cnt gt 0 then bl[nl] = -huge
bu = bnd[1,*]
nl = where(finite(bu) eq 0, cnt)  ; Take non-finite values as 'No bound'
if cnt gt 0 then bu[nl] = huge

key =  n_elements(istate) eq n+1 ? 1 : 0

; If any of A and B is double, we should make the computation in double
if keyword_set(double) then type = size(1d, /TYP) \$
else type = size(A, /TYP) > size(B, /TYP)

;  Initialize flags, etc.
mm = m < n
jj = -1
ifrom5 = 0
iact = 0

; Declaration of variables
x = make_array(n, TYPE=type)               ; solution to be returned
w = make_array(n, TYPE=type)

; The QR solver, la_least_squares, require A to be transposed
; If we favor less memory, we shall to the transpose when calling QR
; but for efficiency, we prefer to transpose A for once:
atr = transpose(a)

; No bl should be greater than huge nor bu smaller than -huge
j = where(bl ge huge, cnt)
if cnt gt 0 then begin
status = 6
if not keyword_set(quiet) then message, /INFO, 'Some lower bounds are excessively large'
return
endif
j = where(bu le -huge, cnt)
if cnt gt 0 then begin
status = 7
if not keyword_set(quiet) then message, /INFO, 'Some upper bounds are excessively small'
return
endif

;  Check consistency of given bounds  bl, bu.
j = where(bl gt -huge and bu lt huge, cnt)
if cnt gt 0 then begin
bdiff = max(bu[j]-bl[j], MIN=bdifmi)
if bdifmi lt 0 then begin
status = 8
if not keyword_set(quiet) then  message, /INFO, 'Inconsistent bounds'
return
endif
if bdiff eq 0 then begin
status = 9
if not keyword_set(quiet) then  message, /INFO, 'No free variables .. check input bounds.'
return
endif
endif

;  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) contains the number of bound variables.  The absolute values of
;  the first nbound 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 eq 0 then begin
nbound = n
nact = 0
istate = -lindgen(nbound+1) - 1
j = where(bl le -huge, cnt)
if cnt gt 0 then istate[j] *= -1    ; no lower bounds
j = where(istate[0:n-1] gt 0 and bu ge huge, cnt, COMPLEM=isbound)
if cnt gt 0 then begin
nbound -= cnt
j = istate[j]
if cnt lt n then istate[0:nbound-1] = istate[isbound]
istate[nbound:n-1] = j           ; unconstrained
istate[n] = nbound
endif
endif else nbound = istate[n]
istate[n] = 0

nact = n - nbound
if nact gt mm then begin
status = 10
if not keyword_set(quiet) then  message, /INFO, 'Too many active variables in starting solution!'
return
endif

; no loop but longer and not clearer (the perforance incidence is tiny)
;j1 = where(istate lt 0, cnt, COMPLEM=j2)
;if cnt gt 0 then begin
;   j = istate[j1]
;   x[j] = bl[j]
;endif
;if cnt lt n then begin
;   j = istate[j2]
;   x[j] = bu[j]
;endif
for k=0, nbound-1 do begin
j = ABS(istate[k]) - 1
if istate[k] lt 0 then x[j] = bl[j] else x[j] = bu[j]
endfor

;  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.
if nbound lt n then begin
kk = istate[nbound:n-1] - 1
j = where(bu[kk] lt huge and bl[kk] gt -huge, cnt, COMPLEM=k)
if cnt gt 0 then x[kk[j]] = (bu[kk[j]]+bl[kk[j]]) / 2 \$
else if cnt lt n-nbound then begin
kk = kk[k]
j1 = where(bu[kk] lt huge, cnt2, COMPLEM=j2)
if cnt2 gt 0 then x[kk[j1]] = bu[kk[j1]] \$
else if cnt2 lt cnt then begin
kk = kk[j2]
j3 = where(bl[kk] gt -huge, cnt3, COMPLEM=j4)
if cnt3 gt 0 then x[kk[j3]] = bl[kk[j3]] \$
else if cnt3 lt cnt2 then x[kk[j4]] = 0
endif
endif
endif

key = nbound lt n ? 1 : 0          ; In cass some variables are unconstrained
;
;!  Compute bnorm, the norm of the data vector b, for reference.
bsq = total(b^2)
bnorm = sqrt(bsq)

;-----------------------------Main Loop---------------------------------
;  Initialization complete.  Begin major loop (Loop A).

for loopa=1,3*n do begin
;  Step 2.
;  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(*,*).

; no gain to make this following computation
;if loopa eq 1 then rsid = b - a # x  \$ ; residuals
;else if nbound lt n then begin
;    j = istate[nbound:n-1] - 1
;    rsid -= a[*,j] # x[j]
;endif
;rsid1 = rsid ; for test

rsid = b - a # x       ; residuals
;rsid2 = rsid ; for test
obj = total(rsid^2)
w = a ## rsid

;  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) le bnorm*eps OR (loopa gt 1 AND nbound eq 0) THEN begin
istate[n] = nbound
w = SQRT(obj)
rnorm = w
return
endif

;  Add the contribution of the active components back into the residual.
if nbound lt n then begin
j = istate[nbound:n-1] - 1
rsid += a[*,j] # x[j]
endif
;rsid3 = rsid ; for test

;  The first iteration in a warm start requires immediate QR.
IF loopa eq 1 AND key ne 0 then GOTO, jump150
;
;  Steps 3, 4.
;  Find the bound element that most wants to be active.
jump120:

worst = 0.
if nbound gt 0 then begin
ks = ABS(istate[0:nbound-1]) - 1
neg = where(istate[0:nbound-1] lt 0, cnt)
if cnt gt 0 then bad[neg] *= -1
iact = ks[it]
endif

;  Test whether the Kuhn-Tucker condition is met.
if worst ge 0d then begin
istate[n] = nbound
w = sqrt(obj)
rnorm = w
return
endif

;  the component  (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 eq jj THEN begin
w[jj] = 0d
GOTO, jump120
ENDIF

;  Step 5.
;  Undo the effect of the new (potentially) active variable on the
;   residual vector.
IF istate[it] gt 0 then bound = bu[iact]
IF istate[it] lt 0 then bound = bl[iact]
rsid += bound * a[*,iact]
;rsid4 = rsid ; for test

;  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-1]
nbound -= 1
nact += 1
istate[nbound] = iact + 1

if mm lt nact then begin
status = 11
if not keyword_set(quiet) then  message, /INFO, ' Too many free variables!'
return
endif

jump150:
;  Step 6.
;  Select 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.

if nact-n+nbound gt 0 then message, 'care ... not tested'
if nact-n+nbound gt 0 then \$
jcl = [jcl[0:nact-n+nbound-1],reverse(istate[nbound:n-1]) - 1] \$
else jcl = reverse(istate[nbound:n-1]) - 1

resq = m ge n ? 0 : -1

zz = la_least_squares(atr[jcl,*], rsid, DOUBLE=double)
if stts eq 0 then stts = check_math(MASK=32) ; dismiss possible underflows

;  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.
;   test.
if resq lt 0d OR \$
(ifrom5 gt 0 AND bu[iact] lt huge AND zz[nact-1] gt bu[iact]) OR  \$
(ifrom5 lt 0 AND bl[iact] gt -huge AND zz[nact-1] lt bl[iact]) THEN begin
;     message, /INFO, 'Warning, this feature has not been tested yet';  ...apparently OK
nbound ++
if  x[iact]-bu[iact] lt 0 then istate[nbound-1] *= -1
nact --
rsid -=  x[iact] * a[*,iact]
ifrom5 = 0
w[iact] = 0d
GOTO, jump120
endif

;  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 ne 0 then jj = -1
ifrom5 = 0

;   Step 7.  Check for strict feasibility of the new QR solution.
j = istate[nbound:nbound+nact-1] - 1
rzz = reverse(zz[0:nact-1])
jmp = where((rzz lt bl[j] and  bl[j] gt -huge) or \$
(rzz gt bu[j] and  bu[j] lt huge), jump)
if jump gt 0 then k1 = jmp + 1

if jump eq 0 then begin ;  New iterate is feasible; back to the top.
x[j] = rzz

endif else begin

;  Steps 8, 9.
alpha = 2d
alf = alpha
for k=k1,nact do begin
j = istate[k+nbound-1] - 1
IF zz[nact-k] gt bu[j] and bu[j] lt huge then alf = (bu[j]-x[j]) / ( zz[nact-k]-x[j])
IF zz[nact-k] lt bl[j] and bl[j] gt -huge  then alf = (bl[j]-x[j]) / ( zz[nact-k]-x[j])
IF alf lt alpha THEN begin
alpha = alf
jj = j
sj = zz[nact-k]-bl[j] ge 0 ? 1d : -1d
endif
endfor

;  Step 10
j = istate[nbound:nbound+nact-1] - 1
x[j] += alpha * (reverse(zz[0:nact-1])-x[j])

;  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)<bl(jj) ).
;   If any other component of  x  is infeasible at this stage, it must
;   be due to roundoff.  Bind every infeasible component and every
;   component at a bound to the appropriate bound.  Correct the
;   residual vector for any variables moved to bounds.  Since at least
;   one variable is removed from the active set in this step, Loop B
;   (Steps 6-11) terminates after at most  nact  steps.
noldb = nbound
for k=1,nact do begin
j = istate[k+noldb-1] - 1
IF (bu[j] lt huge and x[j] gt bu[j]) OR \$
(j eq jj AND sj gt 0d) THEN begin ;  Move x(j) to its upper bound.
x[j] = bu[j]
istate[k+noldb-1] = istate[nbound]
istate[nbound] = j+1
nbound ++
rsid -= bu[j] * a[*,j]
endif else if (bl[j] gt -huge and x[j] le bl[j]) OR \$
(j eq jj AND sj lt 0d) THEN begin ;  Move x(j) to its lower bound.
x[j] = bl[j]
istate[k+noldb-1] = istate[nbound]
istate[nbound] = -j-1
nbound ++
rsid -= bl[j] * a[*,j]
endif
endfor
nact = n - nbound

;  If there are still active variables left repeat the QR; if not,
;    go back to the top.
IF nact gt 0 then GOTO, jump150
endelse ; jump=1
;
endfor

status = 12
if not keyword_set(quiet) then  message, /INFO, 'Failed to converge! '

end ;  bvls_ps

;
PRO Test_BVLS
ncases = 1000
;  ncases = 5
ncols = 10
;  ncols = 3

; Generate artificial data satisfying:  Y = A.beta + noise
seed = 1
beta = randomu(seed, ncols)
beta = 4.0 * (beta - 0.5)
a = randomu(seed, ncases,ncols)
y = dblarr(ncases)
for cas=0l,ncases-1 do begin
e = randomu(seed)
y[cas] = total(a[cas,*] * beta); + 0.1*(e - 0.5)
endfor
huge = (machar()).xmax
undefine, istate
bnd = dblarr(2,ncols)
bnd[0,*] += -1; + (machar()).xmax ;- !VALUES.D_INFINITY
bnd[1,*] += !VALUES.D_INFINITY

;print, 'y=', y
;print, 'a=', a
;print, 'y/a=', y/a

help, a, y
bvls_ps, a, y, bnd, x, istate, LOOP=loopa
bvls_ps, a, y, bnd, x, istate, LOOP=loopa

print, ' Column   Original beta   Solution'
for j=0,ncols-1 do \$
print, FORM='(i5, 2f14.3)', j+1, beta[j], x[j]
print, FORM='(a, i4)', ' No. of iterations = ', loopa

return
END ; PROGRAM Test_BVLS
```
Part of ULySS package