MATLAB: Loosely Constrained Least Squares Solution

least squares solutionlsqlinMATLABnon-linearoptimizationOptimization Toolbox

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 4×21 matrix and d is a 4×1 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

Best Answer

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.

Related Solutions

MATLAB: Are non-linear equality constraints not fulfilled for a constrained optimization problem using “fmincon” and “GlobalSearch”

This behavior is expected, the solver might not be able to satisfy the non-linear constraints completely for all use cases. If the difference is a small value from the bound, the 'TolCon' property can be lowered by doing the following when creating the optimization options:

>> opts = optimoptions(@fmincon,'TolCon',1e-8);

Alternatively you can also try using the alternative algorithms such as ('sqp', 'active-set', 'interior-point') for "fmincon" to see if the constraints are satisfied. More details are available in the following link:

https://www.mathworks.com/help/optim/ug/constrained-nonlinear-optimization-algorithms.html

MATLAB: Quadprog and fmincon provide very different results on the same problem

Why are you doing this? Really, why are you trying to solve it using that kluge of approaches? Go simple. Don't over-complicate things.

n = 4;
AA = [12, -1, 0, 0; ... 
      0, 0, 111, -1];
bb = [0; 0];
x = [1; 10; 1; 61];

Lets see. Minimize norm(x-y), subject to AA*y = bb. Trivial, using the proper tool. In the optimization TB, that is lsqlin. As I show below, my own LSE also gives it directly, found on the file exchange.

Ylsqlin = lsqlin(eye(4),x,[],[],AA,bb)
Minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
Ylsqlin =
         0.834482758578957
          10.0137931029475
         0.549586106124114
          61.0040577797767

Verification:

AA*Ylsqlin - bb
ans =
     -1.77635683940025e-15
                         0

The sum of squares of the difference was:

norm(x - Ylsqlin)^2
ans =
         0.230475348269706

Using my own LSE for comparison, found on the file exchange.

lse(eye(4),x,AA,bb)
ans =
          0.83448275862069
          10.0137931034483
         0.549586106151599
          61.0040577828275

Now, how would it be solved using quadprog? The sum of squares of differences is:

(x-y)'*eye(4)*(x-y)

Expanding, we get

x'*x + y'*y - 2*x'*y

x'*x is a constant, so irrelevant to the optimization, BUT it explains why your objectives are so different.

x'*x
ans =
        3823

To convert this objective to a quadprog usable form, we want it as:

FVAL = 0.5*X'*H*X + f'*X.

H = 2*eye(4);
F = -2*x';

Now solve.

Yqp = quadprog(H,F,[],[],AA,bb)
Minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
Yqp =
         0.834482758599826
          10.0137931031979
         0.549586106137858
          61.0040577813022

It should be no surprise this is the same as lsqlin, but it took more thought to write. Use the proper tools. And now, you want to use fmincon.

fun = @(y) norm(x-y);
x0 = rand(4,1);
Yfmincon = fmincon(fun,x0,[],[],AA,bb)
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
Yfmincon =
         0.834482752432331
           10.013793029188
         0.549586102080044
          61.0040573308849

So esentially identical. And very easy to write in all cases, though lsqlin was by far the simplest.

Anyway, there was absolutely NO need to provide a gradient function to fmincon. (In fact, there is almost never a need to provide a gradient.) A waste of time certainly in this case. Nor even an m-file cost function. The other differences in your results I would put down to simple algebra mistakes that are hardly worth debugging, because you massively over-complicated the problem.

Best Answer

Related Solutions

MATLAB: Are non-linear equality constraints not fulfilled for a constrained optimization problem using “fmincon” and “GlobalSearch”

MATLAB: Quadprog and fmincon provide very different results on the same problem

Related Question