Loosely Constrained Least Squares Solution

Thanks in advance!

I am attempting to find the least squares solution to the matrix equation Ax=b. A, n x m, is a thin matrix, where n>>m, leading to an overdetermined system. Additionally, the system is constrained loosely by the equation Cx=d, where C is a 4x21 matrix and d is a 4x1 vector. Finally, the solutions should fall within the range [-1 1].

For context, I am attempting to optimize a set of Chebyshev polynomial coefficients (hence the solution limit of [-1 1]) defined over the interval [-1 1], expanded to the 20th order, where the polynomial and polynomial first derivative endpoints are pinned to known values (Cx=d). A (n x m) here represents m nominal Chebyshev terms of the first kind (1, x, 2x^2-1, etc.) evaluated at n points between -1 and 1. b represents the value of the original polynomial evaluated at the same n points between -1 and 1, leaving the solution x to be the optimized Chebyshev coefficients.

As such, I applied the following options to lsqlin:

options = optimset('MaxIterations', 1000000, 'TolFun', 1e-5, 'TolCon', 1-e3)

And called lsqlin for each set of coefficients to optimize:

x = lsqlin(C,d,A,b,[],[],-1,1,[],options)

As previously stated, the constraint tolerance is loose, and error on the range of e-3 is acceptable. However, this method of calling lsqlin does not seem to allow for the constraints to have any tolerance. The first derivative endpoints of the resulting polynomials are pinned to the constraint values as desired, but to numerical precision, without any error allowance. This results in lsqlin being unable to find a solution to the problem within the allowed iterations (or even test runs two orders of magnitude higher), and unreasonable results output where coefficients are found to be on the order of e3, not within [-1 1].

Is there something I am missing with the application of constraint tolerances here? Or, is there a way to place a hierarcy on constraints, such that one will be met (coefficients within [-1 1]) at the expense of the other constraints?

Thanks again!

-CH

Answer 1

[John D'Errico]

Your concept of what a constraint tolerance means is wrong. and loose equality constraints are not something you can set.

Instead, you want to set inequality constraints. That is, if you are willing to have these loose constraints, do this:

LooseTol = 1e-3;
CC = [C;-C];
dd = [d + LooseTol ; -d + LooseTol];
lb = -ones(1,size(A,2));
ub = ones(1,size(A,2));

x = lsqlin(A,b,CC,dd,[],[],lb,ub);

Effectively, I have constrained

C*x <= d + LooseTol

C*x >= d - LooseTol

By use of general linear inequality constraints.

[Cy Haukdal]

Thanks a ton for your quick replies, your method certainly seems to be on the right track. I reran my solver with your reccomendations, and there now is some non-zero error for the constrained points. However, it is much smaller than the tolerance I specified, and the [-1 1] solution bounds are still being violated, albeit by a factor of 5 and not 5000.

I specified "looseTol" as 1e-3, and the resulting error was between e-5 and e-12 depending on the run. The solver reached iteration cutoff in all cases, even when resulting error was small:

Maximum number of iterations exceeded; increase options.MaxIter.
To continue solving the problem with the current solution as the
starting point, set x0 = x before calling lsqlin.

I dont see there being a non-linear relationship between constraint tolerance and resulting error at the point constrained. Is that an incorrect assumption?

I will test with some increased error tolerances, but the solver is not fast.

[John D'Errico]

You did not read my comment on Bruno's answer!!!!!

When you set bound constraints on a problem with a VECTOR answer, you NEED TO SPECIFY A VECTOR OF BOUNDS. If you give only scalar bounds, they apply ONLY to the first element of the vector. (And yes, I dislike that behavior. These tools should expand that vector of bounds. That would prevent many errors created by new users.)

A = randn(100,10);
b = randn(100,1);
x = lsqlin(A,b,[],[],[],[],0,1)

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>
x =
        0.0433000144886259
        0.0204015412965199
       0.00845136300344828
        -0.101953066639194
         0.168161930294769
        -0.244447212732897
        -0.122549540051233
       -0.0184344700375564
        0.0799817046345465
         0.131668544182235

So ONLY the very first element of x is constrained.

If you want to constrain ALL of them, then you need to provide a VECTOR of bounds, of the same length as your solution vector.

lsqlin(A,b,[],[],[],[],zeros(10,1),ones(10,1))

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>
ans =
        0.0693449025168689
        0.0306206133110922
      1.82042220055382e-10
      1.83861904547078e-13
         0.203607265647668
      1.69201856659327e-14
      5.22202518593083e-14
      9.46050784216102e-14
         0.059214327404533
         0.104995252144129

And now ALL of the unknowns are bound constrained.

That lsqlin did not satisfy your bound constraints is for the very same reason. All but the very first of your variables had no bound constraints at all applied to them.

[Cy Haukdal]

Sorry, I should have mentioned this explicitly in my comment, I did modify the bounds as you described, such that they are a vector of length m (decreased from 21 to 11 for a lower fidelity test):

lbSoln = chebLim(1)*ones(size(aMat,2),1)

lbSoln =

  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000
  -1.050000000000000

And similarly but positive for upper bound. Bounds are set to +-1.05 since original coefficients were generated poorly and some exceed the [-1 1] bounds slightly, but I dont care about that optimization.

I also fed lsqlin an x0, specifically the nominal coefficients. The solver reached the limits for C1 and C2 here.

disp(coeffArr(1,:,2)')
  -0.394175365020359
   0.459658841701880
   0.092990069979555
   0.020721890270773
   0.004485590022356
   0.000815350813909
   0.000103825497447
  -0.000007832746539
  -0.000011883593189
  -0.000006181149547
  -0.000002077088380

Output:

disp(coeffArrFin(1,:,2)')
  -1.050000000000000
   1.050000000000000
   0.609572679175627
  -0.839229830336653
   0.524146265028548
   0.157789098637682
  -0.487162352832028
   0.185511817332931
   0.099073941787632
  -0.071569577812975
   0.008813392638571

This only happens for a subset of cases, ~50%. In others such as below, optimization is acheived. Note here the minor differences between input and output, this is desired, as the initial conditions are "correct", but require slightly different evaluations at -1 and 1.

Input:

disp(coeffArr(1,:,1)')
  -0.910427905971509
   0.099381108595267
   0.036629347409470
  -0.004698046699548
   0.000234041887144
   0.001169459437181
  -0.001017688304413
   0.000578121566019
  -0.000228213052694
  -0.000064522233648
   0.000104171845374

Output:

disp(coeffArrFin(1,:,1)')
  -0.910427906048462
   0.099381108729015
   0.036629347324175
  -0.004698046664997
   0.000234041885464
   0.001169459427334
  -0.001017688295821
   0.000578121561808
  -0.000228213051369
  -0.000064522233901
   0.000104171845397

Finally, all of the conditions applied via Cx=d are met, in both successful and unsuccessful cases.

[John D'Errico]

I don't understand. Before you were saying the bounds were not being satisfied. Clearly they were in the cases you show.

In those examples, we see one case where the solution you expected was not the one that was found. In another case, the solution is the same.

What we do not know is if the linear system may be well posed or not. If it is rank deficient, then you may not get a unique solution. In fact, you can get infinitely many solutions, so the solution you expect need not be the one you get.

It is a little hard to know what is the issue here. I would do things like test the condition number of the matrix A. Is it large? (I.e., near 1/eps?) If so, then problems will exist. Is the solution you think to be "correct" really a valid solution? Are the constraints active at that set? What would the unconstrained problem yield as a solution?

[Cy Haukdal]

To be more specific, a Chebyshev polynomial should both be characterized by coefficients and values within [-1 1].

In this case, the bounds applied via Cx=d are met, the bounds on the soultion values are met [-1 1], but the problem is incorrectly solved, indicated both by relation to expected solution and values outside of [-1 1].

To answer your questions:

• The problem is not well posed if multiple solutions exist, although since lsqlin is exiting based on iterations, its possible the answers are not solutions
• The matrix A is not rank deficient
• The condition number of A is 2.9
• It is unknown if the solutions are valid, but it is assumed, since there are successful cases
• "correct" solutions meet all constraints imposed on Ax=b
• Unconstrained problem A\b produces the correct/expected solution in all cases, but with end conditions which are insufficient for my purposes. Also implies solutions are valid.

Even in the above runs, the error at constrained points is orders of magnitude below the tolerance. I think the issue still lies here.

Answer 2

[Bruno Luong]

alpha = some positive value; % larget alpha, norm(C*x-d) gets smaller

E = [A; alpha*C];
f = [b; alpha*d];

x = lsqlin(E,f,[],[],[],[],-1,1,[],options)

[John D'Errico]

An important point was missed here. When x is a vector, but the bounds are supplied only as a scalar, then the bounds only apply to the FIRST element of the vector. (I never liked that, but it is true. Far too many people make this mistake.)

So this is not valid for vector x.

x = lsqlin(E,f,[],[],[],[],-1,1,[],options)

Instead, use:

n = size(A,2);
lb = -ones(n,1);
ub = ones(n,1);
x = lsqlin(E,f,[],[],[],[],lb,ub,[],options);

If a scalar only is supplied, then consider this simple example:

lsqlin(randn(100,10),randn(100,1),[],[],[],[],0,1)

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>
ans =
       4.4012355548152e-10
        0.0377174311402519
       -0.0981597231592627
        -0.102733917267231
        0.0117681433517298
       -0.0666227293207337
         0.127194321772405
         -0.12742206265685
        -0.101845196324251
        0.0025344777053746

As you can see, only the first element was constrained to have a lower bound of 0. The others were fully unconstrained.

[Bruno Luong]

You might be right John. I just copy OP's bounds and focussing on the change of cost function and equality constraints.

[John D'Errico]

Yes. I realize that you surely know the difference. This was more for the benefit of the OP, or others who would see your response.