Solve constrained linear least-squares problems
Linear least-squares solver with bounds or linear constraints.
Solves least-squares curve fitting problems of the form
Note
lsqlin
applies only to the solver-based approach. For a discussion
of the two optimization approaches, see First Choose Problem-Based or Solver-Based Approach.
finds
the minimum for x
= lsqlin(problem
)problem
, a structure described in problem
. Create the problem
structure using dot notation or
the struct
function. Or create a
problem
structure from an OptimizationProblem
object
by using prob2struct
.
[
, for any input arguments described
above, returns:x
,resnorm
,residual
,exitflag
,output
,lambda
]
= lsqlin(___)
The squared 2-norm of the residual resnorm =
The residual residual = C*x - d
A value exitflag
describing the
exit condition
A structure output
containing information
about the optimization process
A structure lambda
containing the
Lagrange multipliers
The factor ½ in the definition of the problem affects the
values in the lambda
structure.
Find the x
that minimizes the norm of C*x - d
for an overdetermined problem with linear inequality constraints.
Specify the problem and constraints.
C = [0.9501 0.7620 0.6153 0.4057 0.2311 0.4564 0.7919 0.9354 0.6068 0.0185 0.9218 0.9169 0.4859 0.8214 0.7382 0.4102 0.8912 0.4447 0.1762 0.8936]; d = [0.0578 0.3528 0.8131 0.0098 0.1388]; A = [0.2027 0.2721 0.7467 0.4659 0.1987 0.1988 0.4450 0.4186 0.6037 0.0152 0.9318 0.8462]; b = [0.5251 0.2026 0.6721];
Call lsqlin
to solve the problem.
x = lsqlin(C,d,A,b)
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.
x = 4×1
0.1299
-0.5757
0.4251
0.2438
Find the x
that minimizes the norm of C*x - d
for an overdetermined problem with linear equality and inequality constraints and bounds.
Specify the problem and constraints.
C = [0.9501 0.7620 0.6153 0.4057 0.2311 0.4564 0.7919 0.9354 0.6068 0.0185 0.9218 0.9169 0.4859 0.8214 0.7382 0.4102 0.8912 0.4447 0.1762 0.8936]; d = [0.0578 0.3528 0.8131 0.0098 0.1388]; A =[0.2027 0.2721 0.7467 0.4659 0.1987 0.1988 0.4450 0.4186 0.6037 0.0152 0.9318 0.8462]; b =[0.5251 0.2026 0.6721]; Aeq = [3 5 7 9]; beq = 4; lb = -0.1*ones(4,1); ub = 2*ones(4,1);
Call lsqlin
to solve the problem.
x = lsqlin(C,d,A,b,Aeq,beq,lb,ub)
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.
x = 4×1
-0.1000
-0.1000
0.1599
0.4090
This example shows how to use nondefault options for linear least squares.
Set options to use the 'interior-point'
algorithm and to give iterative display.
options = optimoptions('lsqlin','Algorithm','interior-point','Display','iter');
Set up a linear least-squares problem.
C = [0.9501 0.7620 0.6153 0.4057 0.2311 0.4564 0.7919 0.9354 0.6068 0.0185 0.9218 0.9169 0.4859 0.8214 0.7382 0.4102 0.8912 0.4447 0.1762 0.8936]; d = [0.0578 0.3528 0.8131 0.0098 0.1388]; A = [0.2027 0.2721 0.7467 0.4659 0.1987 0.1988 0.4450 0.4186 0.6037 0.0152 0.9318 0.8462]; b = [0.5251 0.2026 0.6721];
Run the problem.
x = lsqlin(C,d,A,b,[],[],[],[],[],options)
Iter Fval Primal Infeas Dual Infeas Complementarity 0 -7.687420e-02 1.600492e+00 6.150431e-01 1.000000e+00 1 -7.687419e-02 8.002458e-04 3.075216e-04 2.430833e-01 2 -3.162837e-01 4.001229e-07 1.537608e-07 5.945636e-02 3 -3.760545e-01 2.000616e-10 2.036997e-08 1.370933e-02 4 -3.912129e-01 1.000866e-13 1.006816e-08 2.548273e-03 5 -3.948062e-01 2.220446e-16 2.955102e-09 4.295807e-04 6 -3.953277e-01 2.775558e-17 1.237758e-09 3.102850e-05 7 -3.953581e-01 2.775558e-17 1.645862e-10 1.138719e-07 8 -3.953582e-01 2.775558e-17 2.399608e-13 5.693290e-11 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.
x = 4×1
0.1299
-0.5757
0.4251
0.2438
Obtain and interpret all lsqlin
outputs.
Define a problem with linear inequality constraints and bounds. The problem is overdetermined because there are four columns in the C
matrix but five rows. This means the problem has four unknowns and five conditions, even before including the linear constraints and bounds.
C = [0.9501 0.7620 0.6153 0.4057 0.2311 0.4564 0.7919 0.9354 0.6068 0.0185 0.9218 0.9169 0.4859 0.8214 0.7382 0.4102 0.8912 0.4447 0.1762 0.8936]; d = [0.0578 0.3528 0.8131 0.0098 0.1388]; A = [0.2027 0.2721 0.7467 0.4659 0.1987 0.1988 0.4450 0.4186 0.6037 0.0152 0.9318 0.8462]; b = [0.5251 0.2026 0.6721]; lb = -0.1*ones(4,1); ub = 2*ones(4,1);
Set options to use the 'interior-point'
algorithm.
options = optimoptions('lsqlin','Algorithm','interior-point');
The 'interior-point'
algorithm does not use an initial point, so set x0
to []
.
x0 = [];
Call lsqlin
with all outputs.
[x,resnorm,residual,exitflag,output,lambda] = ...
lsqlin(C,d,A,b,[],[],lb,ub,x0,options)
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.
x = 4×1
-0.1000
-0.1000
0.2152
0.3502
resnorm = 0.1672
residual = 5×1
0.0455
0.0764
-0.3562
0.1620
0.0784
exitflag = 1
output = struct with fields:
message: '...'
algorithm: 'interior-point'
firstorderopt: 4.3374e-11
constrviolation: 0
iterations: 6
linearsolver: 'dense'
cgiterations: []
lambda = struct with fields:
ineqlin: [3x1 double]
eqlin: [0x1 double]
lower: [4x1 double]
upper: [4x1 double]
Examine the nonzero Lagrange multiplier fields in more detail. First examine the Lagrange multipliers for the linear inequality constraint.
lambda.ineqlin
ans = 3×1
0.0000
0.2392
0.0000
Lagrange multipliers are nonzero exactly when the solution is on the corresponding constraint boundary. In other words, Lagrange multipliers are nonzero when the corresponding constraint is active. lambda.ineqlin(2)
is nonzero. This means that the second element in A*x
should equal the second element in b
, because the constraint is active.
[A(2,:)*x,b(2)]
ans = 1×2
0.2026 0.2026
Now examine the Lagrange multipliers for the lower and upper bound constraints.
lambda.lower
ans = 4×1
0.0409
0.2784
0.0000
0.0000
lambda.upper
ans = 4×1
0
0
0
0
The first two elements of lambda.lower
are nonzero. You see that x(1)
and x(2)
are at their lower bounds, -0.1
. All elements of lambda.upper
are essentially zero, and you see that all components of x
are less than their upper bound, 2
.
C
— Multiplier matrixMultiplier matrix, specified as a matrix of doubles. C
represents
the multiplier of the solution x
in the expression C*x
- d
. C
is M
-by-N
,
where M
is the number of equations, and N
is
the number of elements of x
.
Example: C = [1,4;2,5;7,8]
Data Types: double
d
— Constant vectorConstant vector, specified as a vector of doubles. d
represents
the additive constant term in the expression C*x - d
. d
is M
-by-1
,
where M
is the number of equations.
Example: d = [5;0;-12]
Data Types: double
A
— Linear inequality constraintsLinear inequality constraints, specified as a real matrix. A
is
an M
-by-N
matrix, where M
is
the number of inequalities, and N
is the number
of variables (number of elements in x0
). For
large problems, pass A
as a sparse matrix.
A
encodes the M
linear
inequalities
A*x <= b
,
where x
is the column vector of N
variables x(:)
,
and b
is a column vector with M
elements.
For example, to specify
x1 + 2x2 ≤
10
3x1 +
4x2 ≤ 20
5x1 +
6x2 ≤ 30,
enter these constraints:
A = [1,2;3,4;5,6]; b = [10;20;30];
Example: To specify that the x components sum to 1 or less, use A =
ones(1,N)
and b = 1
.
Data Types: double
b
— Linear inequality constraintsLinear inequality constraints, specified as a real vector. b
is
an M
-element vector related to the A
matrix.
If you pass b
as a row vector, solvers internally
convert b
to the column vector b(:)
.
For large problems, pass b
as a sparse vector.
b
encodes the M
linear
inequalities
A*x <= b
,
where x
is the column vector of N
variables x(:)
,
and A
is a matrix of size M
-by-N
.
For example, consider these inequalities:
x1
+ 2x2 ≤
10
3x1
+ 4x2 ≤
20
5x1
+ 6x2 ≤
30.
Specify the inequalities by entering the following constraints.
A = [1,2;3,4;5,6]; b = [10;20;30];
Example: To specify that the x components sum to 1 or less, use A =
ones(1,N)
and b = 1
.
Data Types: double
Aeq
— Linear equality constraintsLinear equality constraints, specified as a real matrix. Aeq
is
an Me
-by-N
matrix, where Me
is
the number of equalities, and N
is the number of
variables (number of elements in x0
). For large
problems, pass Aeq
as a sparse matrix.
Aeq
encodes the Me
linear
equalities
Aeq*x = beq
,
where x
is the column vector of N
variables x(:)
,
and beq
is a column vector with Me
elements.
For example, to specify
x1 + 2x2 +
3x3 = 10
2x1 +
4x2 + x3 =
20,
enter these constraints:
Aeq = [1,2,3;2,4,1]; beq = [10;20];
Example: To specify that the x components sum to 1, use Aeq = ones(1,N)
and
beq = 1
.
Data Types: double
beq
— Linear equality constraintsLinear equality constraints, specified as a real vector. beq
is
an Me
-element vector related to the Aeq
matrix.
If you pass beq
as a row vector, solvers internally
convert beq
to the column vector beq(:)
.
For large problems, pass beq
as a sparse vector.
beq
encodes the Me
linear
equalities
Aeq*x = beq
,
where x
is the column vector of N
variables
x(:)
, and Aeq
is a matrix of size
Me
-by-N
.
For example, consider these equalities:
x1
+ 2x2 +
3x3 =
10
2x1
+ 4x2 +
x3 =
20.
Specify the equalities by entering the following constraints.
Aeq = [1,2,3;2,4,1]; beq = [10;20];
Example: To specify that the x components sum to 1, use Aeq = ones(1,N)
and
beq = 1
.
Data Types: double
lb
— Lower bounds[]
(default) | real vector or arrayLower bounds, specified as a vector or array of doubles. lb
represents
the lower bounds elementwise in lb
≤ x
≤ ub
.
Internally, lsqlin
converts an array lb
to
the vector lb(:)
.
Example: lb = [0;-Inf;4]
means x(1)
≥ 0
, x(3) ≥ 4
.
Data Types: double
ub
— Upper bounds[]
(default) | real vector or arrayUpper bounds, specified as a vector or array of doubles. ub
represents
the upper bounds elementwise in lb
≤ x
≤ ub
.
Internally, lsqlin
converts an array ub
to
the vector ub(:)
.
Example: ub = [Inf;4;10]
means x(2)
≤ 4
, x(3) ≤ 10
.
Data Types: double
x0
— Initial point[]
(default) | real vector or arrayInitial point for the solution process, specified as a real vector or array. The
'trust-region-reflective'
and 'active-set'
algorithms use x0
(optional).
If you do not specify x0
for the
'trust-region-reflective'
or 'active-set'
algorithm, lsqlin
sets x0
to the zero vector. If
any component of this zero vector x0
violates the bounds,
lsqlin
sets x0
to a point in the interior of
the box defined by the bounds.
Example: x0 = [4;-3]
Data Types: double
options
— Options for lsqlin
optimoptions
| structure such as created by optimset
Options for lsqlin
, specified as the output of the optimoptions
function or as a structure such as created by
optimset
.
Some options are absent from the
optimoptions
display. These options appear in italics in the following
table. For details, see View Options.
All Algorithms
| Choose the algorithm:
The The
When the problem has no constraints,
If you have a large number of linear constraints and
not a large number of variables, try the For more information on choosing the algorithm, see Choosing the Algorithm. |
Diagnostics | Display diagnostic information about the function to
be minimized or solved. The choices are |
Display | Level of display returned to the command line.
The
|
MaxIterations | Maximum number of iterations allowed, a positive integer. The default value is
For
|
trust-region-reflective
Algorithm Options
FunctionTolerance | Termination tolerance on the function value, a positive scalar.
The default is For |
JacobianMultiplyFcn | Jacobian multiply function, specified
as a function handle. For large-scale structured problems, this function
should compute the Jacobian matrix product W = jmfun(Jinfo,Y,flag) where
In each case, See Jacobian Multiply Function with Linear Least Squares for an example. For |
MaxPCGIter | Maximum number of PCG (preconditioned
conjugate gradient) iterations, a positive scalar. The default is |
OptimalityTolerance | Termination tolerance on the first-order optimality, a positive
scalar. The default is For |
PrecondBandWidth | Upper bandwidth of preconditioner
for PCG (preconditioned conjugate gradient). By default, diagonal
preconditioning is used (upper bandwidth of 0). For some problems,
increasing the bandwidth reduces the number of PCG iterations. Setting |
SubproblemAlgorithm | Determines how the iteration step
is calculated. The default, |
TolPCG | Termination tolerance on the PCG
(preconditioned conjugate gradient) iteration, a positive scalar.
The default is |
TypicalX | Typical |
interior-point
Algorithm Options
ConstraintTolerance | Tolerance on the constraint violation, a positive scalar. The default is
For |
LinearSolver | Type of internal linear solver in algorithm:
|
OptimalityTolerance | Termination tolerance on the first-order optimality, a positive
scalar. The default is For |
StepTolerance | Termination tolerance on For
|
'active-set'
Algorithm Options
ConstraintTolerance | Tolerance on the constraint violation, a positive scalar. The default
value is For |
ObjectiveLimit | A tolerance (stopping criterion) that is a
scalar. If the objective function value goes below
|
OptimalityTolerance | Termination tolerance on the first-order
optimality, a positive scalar. The default value is For |
StepTolerance | Termination tolerance on For |
problem
— Optimization problemOptimization problem, specified as a structure with the following fields.
| Matrix multiplier in the term C*x
- d |
| Additive constant in the term C*x
- d |
| Matrix for linear inequality constraints |
| Vector for linear inequality constraints |
| Matrix for linear equality constraints |
| Vector for linear equality constraints |
lb | Vector of lower bounds |
ub | Vector of upper bounds |
| Initial point for x |
| 'lsqlin' |
| Options created with optimoptions |
Data Types: struct
x
— SolutionSolution, returned as a vector that minimizes the norm of C*x-d
subject
to all bounds and linear constraints.
resnorm
— Objective valueObjective value, returned as the scalar value norm(C*x-d)^2
.
residual
— Solution residualsSolution residuals, returned as the vector C*x-d
.
exitflag
— Algorithm stopping conditionAlgorithm stopping condition, returned as an integer identifying
the reason the algorithm stopped. The following lists the values of exitflag
and
the corresponding reasons lsqlin
stopped.
| Change in the residual was smaller than the specified tolerance
|
| Step size smaller than
|
| Function converged to a solution
|
| Number of iterations exceeded |
| The problem is infeasible. Or, for the |
-3 | The problem is unbounded. |
| Ill-conditioning prevents further optimization. |
| Unable to compute a step direction. |
The exit message for the interior-point
algorithm
can give more details on the reason lsqlin
stopped,
such as exceeding a tolerance. See Exit Flags and Exit Messages.
output
— Solution process summarySolution process summary, returned as a structure containing information about the optimization process.
| Number of iterations the solver took. |
| One of these algorithms:
For an unconstrained problem, |
| Constraint violation that is positive
for violated constraints (not returned for the
|
| Exit message. |
| First-order optimality at the solution. See First-Order Optimality Measure. |
linearsolver | Type of internal linear solver,
|
| Number of conjugate gradient iterations
the solver performed. Nonempty only for the |
See Output Structures.
lambda
— Lagrange multipliersLagrange multipliers, returned as a structure with the following fields.
| Lower bounds |
| Upper bounds |
| Linear inequalities |
| Linear equalities |
For problems with no constraints, you can use mldivide
(matrix left division).
When you have no constraints, lsqlin
returns x = C\d
.
Because the problem being solved is always convex, lsqlin
finds
a global, although not necessarily unique, solution.
If your problem has many linear constraints and few variables, try using the
'active-set'
algorithm. See Quadratic Programming with Many Linear Constraints.
Better numerical results are likely if you specify
equalities explicitly, using Aeq
and beq
,
instead of implicitly, using lb
and ub
.
The trust-region-reflective
algorithm
does not allow equal upper and lower bounds. Use another algorithm
for this case.
If the specified input bounds for a problem are inconsistent,
the output x
is x0
and the outputs resnorm
and residual
are []
.
You can solve some large structured problems, including
those where the C
matrix is too large to fit in
memory, using the trust-region-reflective
algorithm
with a Jacobian multiply function. For information, see trust-region-reflective Algorithm Options.
This method is a subspace trust-region method based on the interior-reflective Newton method described in [1]. Each iteration involves the approximate solution of a large linear system using the method of preconditioned conjugate gradients (PCG). See Trust-Region-Reflective Least Squares, and in particular Large Scale Linear Least Squares.
The 'interior-point'
algorithm is based on
the quadprog
'interior-point-convex'
algorithm.
See Linear Least Squares: Interior-Point or Active-Set.
The 'active-set'
algorithm is based on the
quadprog
'active-set'
algorithm. For more information, see Linear Least Squares: Interior-Point or Active-Set and active-set quadprog Algorithm.
[1] Coleman, T. F. and Y. Li. “A Reflective Newton Method for Minimizing a Quadratic Function Subject to Bounds on Some of the Variables,” SIAM Journal on Optimization, Vol. 6, Number 4, pp. 1040–1058, 1996.
[2] Gill, P. E., W. Murray, and M. H. Wright. Practical Optimization, Academic Press, London, UK, 1981.
The Optimize Live Editor task provides a visual interface for lsqlin
.
You have a modified version of this example. Do you want to open this example with your edits?