BVLS and Kuhn-Tucker condition

[svlad]

I'm using an f2c'd version of the BVLS Fortran code to solve a bounded
linear, least-squares problem. I've run across a scenario where the
BVLS code exits via the Kuhn-Tucker condition when the current solution
is far from optimal. I have the following very simple system:

| 75.2193 -2.59184 -2.25085 |
A = | -2.59184 3.46673 -0.0311589 |
| -2.25085 -0.0311589 2.83563 |

b = [ 94.4536 0.594381 1.41701 ]

with all variables bounded to a range [0, 1]

For a "cold-start" with all variables set to their lower-bound, in the
first iteration the negative gradient vector is:

w = [7100.01 -242.793 -208.601 ]

Transfering the sign of the bound (-1 for lower-bound, +1 for
upper-bound), the first variable has the most negative value and
becomes the best candidate for activation. The QR solution for this
variable is 1.25227. The solution is infeasible resulting in the
variable getting clamped to its upper-bound: 1.0. There are no more
active variables so the algorithm begins the second iteration.

In the second iteration the negative gradient vector is now:

w = [1430.28 -38.9209 -32.9922 ]

The first variable is at its upper-bound and the last two are at their
lower-bound. After transfering the sign of the bounds (-1 for
lower-bound, +1 for upper-bound) we are left with no variables with a
negative gradient value and the algorithm terminates via Kuhn-Tucker
with the following solution:

x = [1.0 0.0 0.0 ]

Solving the same system with a simple Gauss-Seidel iteration clamping
each variable to its bound at each iteration gives the following
result:

x = [ 1.0 0.9280734 1.0]

Clearly the BVLS solution is nowhere near optimal.

My understanding of the Kuhn-Tucker condition and the gradient
computation is not very clear. Intuitively, I would think that the
gradient vector is indicating the direction each variable's value needs
to move in order to reduce the residual. The Kuhn-Tucker condition then
provides a method for determining when all bound variables would have
to move in a direction that takes them beyond their bound in that
direction (e.g. already at bound of 1.0 and wanting to become larger in
value).

An unconstrained solution of the same system yields the following
result:

x = [ 1.3439983 1.1904674 1.5796286 ]

The solution (which has a residual of 0.0) has each variable beyond the
upper-bound of 1.0. I would think the gradient vector for the system
should have three positive values? Is there something inherently wrong
with the gradient vector computation or is this simply a limitation or
shortcoming of this method?

Any insights would be *greatly* appreciated.

Thanks!

[DTB]

There are at least two versions of a code named BVLS in the open
software literature. They do the same job. Who are the authors of the
version you are using? Stark and Parker or Lawson and Hanson?
DTB

[svlad]

I am using the Stark and Parker version.

[Peter Spellucci]

In article <1164070797....@h48g2000cwc.googlegroups.com>,

WHICH IS CORRECT

>Solving the same system with a simple Gauss-Seidel iteration clamping
>each variable to its bound at each iteration gives the following
>result:
>
>x = [ 1.0 0.9280734 1.0]
>
>Clearly the BVLS solution is nowhere near optimal.
>
>My understanding of the Kuhn-Tucker condition and the gradient
>computation is not very clear. Intuitively, I would think that the
>gradient vector is indicating the direction each variable's value needs
>to move in order to reduce the residual. The Kuhn-Tucker condition then
>provides a method for determining when all bound variables would have
>to move in a direction that takes them beyond their bound in that
>direction (e.g. already at bound of 1.0 and wanting to become larger in
>value).
>
>An unconstrained solution of the same system yields the following
>result:
>
>x = [ 1.3439983 1.1904674 1.5796286 ]
>
>The solution (which has a residual of 0.0) has each variable beyond the
>upper-bound of 1.0. I would think the gradient vector for the system
>should have three positive values? Is there something inherently wrong
>with the gradient vector computation or is this simply a limitation or
>shortcoming of this method?
>
>Any insights would be *greatly* appreciated.
>
>Thanks!
>

YOU HAD BETTER DONE TO TEST THE ERROR CRITERION:
HERE IS THE ITERATION WITH YOUR "OPTIMAL" GAUSS SEIDEL INPUT:
1 FX= .5715008D+03 UPSI= .00D+00 B2N= .00D+00 UMI= -.12D+03 NR 3 SI-1
2 FX= .3648776D+03 UPSI= .41D+00 B2N= .99D+03 UMI= .00D+00 NR 4 SI 1
3 FX= .3053753D+03 UPSI= .58D-01 B2N= .17D+04 UMI= .00D+00 NR 4 SI 1
4 FX= .8011123D+02 UPSI= .17D+00 B2N= .40D+03 UMI= .00D+00 NR 5 SI 1
5 FX= .5257758D+02 UPSI= .19D+00 B2N= .23D-12 UMI= -.18D+01 NR 3 SI-1
6 FX= .7585253D+02 UPSI= .17D+00 B2N= .36D+03 UMI= .00D+00 NR 1 SI-1
7 FX= .9157528D+02 UPSI= .44D+00 B2N= .13D+03 UMI= .00D+00 NR 4 SI 1
8 FX= .1407725D+03 UPSI= .33D+00 B2N= .45D-12 UMI= .00D+00 NR 3 SI-1

NORM OF GRAD(F) .2862373D+04
LAGRANGIAN VIOLATION .2009718D-13
FEASIBILITY VIOLATION .1110223D-15
DUAL FEASIBILITY VIOLATION .0000000D+00
OPTIMIZER RUNTIME SEC'S .1300000E+00

OPTIMAL VALUE OF F = .393563497730441D+03

OPTIMAL SOLUTION X =
.100000000000000D+01 .144198888940572D-15 -.111022302462516D-15

(FOR MY CONVENIENCE, I SOLVED IT WITH MY CODE DONLP2)
HTH
PETER

[svlad]

So it appears that I am misunderstanding something in regards to the
measurement of error or "quality" of the solution. In my application, I
am simply trying to align one set of points against a set of target
positions. The points move linearly and are controlled by the variables
for which I am trying to find a set of optimal values minimizing the
distance between the set of points.

It is in that sense that the bounded Gauss-Seidel iterations produce a
more "optimal" fit to the target points. In the sample system I posted,
there are twelve points whose positions are controlled by three
variables. The sum of the distances between the points with the BVLS
solution (1, 0, 0) is 7.25. The sum of the distances between the points
with the Gauss-Seidel bounded iterations is 5.35. I wish I could post
an image because the disparity between the two solutions is visually
apparent.

I'm very confused because I thought that the gradient vector was
measuring which variable would have the greatest impact in reducing the
error in the system i.e. minimizing the squared distances between the
points. What is the gradient vector actually measuring? Any ideas on
how I can solve this in the sense I desire?

Thanks!

[mecej4]

The comments at the top of BVLS.F clearly state that the code minimizes
the 2-norm ||A x - b||_2 while honoring the bound constraints. If that
2-norm is not what you want minimized, your usage of BVLS is incorrect,
and no conclusions about the result should be drawn: there is no basis
for saying that "..the BVLS solution is nowhere near optimal..", since
"optimal" has not been defined in your posts in terms of A, x and b, nor
have you stated how the "squared distances between the points" are
related to A, x and b.

-- mecej4

[Peter Spellucci]

In article <1164147682.7...@k70g2000cwa.googlegroups.com>,

this sounds as if you didn't want to optimize the euclidean distance
of A*x and b subject to x in the unit box, but something quite different.
the gradient points in teh direction of steepest acsent , steepest in the
sense of euclidean geometry!
maybe your tru problem is something like
sum_{i = 1 to N} dist(x(var)(i),xdata(i)) with respect to "var"? (!)
this is quite different from
dist(x(var), xdata)
where x and xdata are now vectors composed from N subvectors.
you can transform problem (!) into an constrained smooth optimization problem
with quadratic constraints:
dist^2( xvar(i),xdata(i) ) <= epsilon(i)
sum_{i=1 to N} epsilon(i) = min
hth
peter

[svlad]

Here is a description of the problem I am trying to solve:

http://groups.google.com/group/sci.math.num-analysis/browse_thread/thread/4b377fb2779cd504

Clearly in the unbounded case the minimization of the gradient norm of
the distance squared function is what I am after. What has been pointed
out to me thus far is that in the face of the bounds this is no longer
true for my application.

Peter Spellucci wrote:
> this sounds as if you didn't want to optimize the euclidean distance
> of A*x and b subject to x in the unit box, but something quite different.
> the gradient points in teh direction of steepest acsent , steepest in the
> sense of euclidean geometry!
> maybe your tru problem is something like
> sum_{i = 1 to N} dist(x(var)(i),xdata(i)) with respect to "var"? (!)
> this is quite different from
> dist(x(var), xdata)
> where x and xdata are now vectors composed from N subvectors.
> you can transform problem (!) into an constrained smooth optimization problem
> with quadratic constraints:
> dist^2( xvar(i),xdata(i) ) <= epsilon(i)
> sum_{i=1 to N} epsilon(i) = min

I've just gone off the edge of my very small map of knowledge in this
domain. Any pointers on how to set the problem up in this manner? Also
will these NLP techniques allow me to implement mutual exclusivity
constraints in addition to the bound constraints (e.g. variable x[1]
and x[5] cannot both be > 0 where x[1] and x[5] must both be [0,1])?

Thanks!

[Ray Koopman]

Check your Gauss-Seidel routine. Mine gives x = {1, 0, 0}.

[svlad]

I did a quick test in SciLab using the following algorithm:

function [x, toler, k] = gauss_seidel_bounded(A, b, bl, bu, x0, eps,
maxIter)

n = size(A, 'c');
x = x0';
k = 0;
toler = 1;
while ((k < maxIter) & (toler > eps))

xold = x;
for i=1:n

sigma = 0.0;
for j=1:n
sigma = sigma + (A(i, j) * x(j));
end

x(i) = (1 / A(i,i)) * (b(i) + (A(i,i) * x(i)) - sigma);

if(x(i) < bl(i))
x(i) = bl(i);
end

if(x(i) > bu(i))
x(i) = bu(i);
end

end

toler = norm((x-xold));
k = k + 1;

end

endfunction

The algorithm terminates after 3 iterations with: x = [1.0 0.9280665
1.0] and a residual of [-23.887331, 0.0, -0.8615518].

Is this algorithm different than what you used?

If I solve the same system unbounded I get: x = [1.3439579 1.190431
1.5797378] with a residual of [0, 0, 0]. This answer seems to line up
with what any old direct solution method would give.

[Peter Spellucci]

In article <1164677865....@h54g2000cwb.googlegroups.com>,

this reveals the misunderstanding:
you minimize
(1/2)*x'*A*x-b'*x subject to x in [bl,bu]
Ray and I (and you with BVLS) minimized

(1/2)||Ax-b||_2^2 = (1/2)*x'*A'*A*x-b'*A*x subject to x in [bl,bu]

this is quite different!
hth
peter