Solves an lsei inverse problem (Least Squares with Equality and Inequality
Constraints)
min(||Ax-b||^2)
subject to
Ex=f
Gx>=h
Uses either subroutine lsei (FORTRAN) from the LINPACK package, or
solve.QP from R-package quadprog.
In case the equality constraints Ex=f cannot be satisfied, a
generalized inverse solution residual vector length is obtained for f-Ex.
This is the minimal length possible for ||f-Ex||^2.
Usage
lsei (A = NULL, B = NULL, E = NULL, F = NULL, G = NULL, H = NULL,
Wx = NULL, Wa = NULL, type = 1, tol = sqrt(.Machine$double.eps),
tolrank = NULL, fulloutput = FALSE, verbose = TRUE)
Arguments
A
numeric matrix containing the coefficients of the quadratic
function to be minimised, ||Ax-B||^2; if the columns of A have
a names attribute, they will be used to label the output.
B
numeric vector containing the right-hand side of the quadratic
function to be minimised.
E
numeric matrix containing the coefficients of the equality
constraints, Ex=F; if the columns of E have a names attribute,
and the columns of A do not, they will be used to label the output.
F
numeric vector containing the right-hand side of the equality
constraints.
G
numeric matrix containing the coefficients of the inequality
constraints, Gx>=H; if the columns of G have a names
attribute, and the columns of A and E do not, they will be
used to label the output.
H
numeric vector containing the right-hand side of the inequality
constraints.
Wx
numeric vector with weighting coefficients of unknowns (length
= number of unknowns).
Wa
numeric vector with weighting coefficients of the quadratic
function (Ax-B) to be minimised (length = number of number of rows of A).
type
integer code determining algorithm to use 1=lsei,
2=solve.QP from R-package quadprog (see note).
tol
tolerance (for singular value decomposition, equality and
inequality constraints).
tolrank
only used if type = 1; if not NULL then
tolrank should be a two-valued vector containing the rank
determination tolerance for the equality constraint equations
(1st value) and for the reduced least squares equations (2nd value).
fulloutput
if TRUE, also returns the covariance matrix of
the solution and the rank of the equality constraints - only used
if type = 1.
verbose
logical to print error messages.
Value
a list containing:
X
vector containing the solution of the least squares problem.
residualNorm
scalar, the sum of absolute values of residuals of
equalities and violated inequalities.
solutionNorm
scalar, the value of the minimised quadratic function
at the solution, i.e. the value of ||Ax-b||^2.
IsError
logical, TRUE if an error occurred.
type
the string "lsei", such that how the solution was obtained
can be traced.
covar
covariance matrix of the solution; only returned if
fulloutput = TRUE.
RankEq
rank of the equality constraint matrix.; only returned if
fulloutput = TRUE.
RankApp
rank of the reduced least squares problem (approximate
equations); only returned if fulloutput = TRUE.
Note
See comments in the original code for more details; these comments are
included in the ‘docs’ subroutine of the package.
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
References
K. H. Haskell and R. J. Hanson, An algorithm for linear least squares problems
with equality and nonnegativity constraints, Report SAND77-0552,
Sandia Laboratories, June 1978.
K. H. Haskell and R. J. Hanson, Selected algorithms for the linearly
constrained least squares problem - a users guide,
Report SAND78-1290, Sandia Laboratories,August 1979.
K. H. Haskell and R. J. Hanson, An algorithm for linear least squares
problems with equality and nonnegativity constraints,
Mathematical Programming 21 (1981), pp. 98-118.
R. J. Hanson and K. H. Haskell, Two algorithms for the linearly constrained
least squares problem, ACM Transactions on Mathematical Software,
September 1982.
Berwin A. Turlach R and Andreas Weingessel (2007). quadprog: Functions to
solve Quadratic Programming Problems. R package version 1.4-11. S original
by Berwin A. Turlach R port by Andreas Weingessel.
See Also
ldei, linp,
solve.QR the original function from package quadprog.
Examples
# ------------------------------------------------------------------------------
# example 1: polynomial fitting
# ------------------------------------------------------------------------------
x <- 1:5
y <- c(9, 8, 6, 7, 5)
plot(x, y, main = "Polynomial fitting, using lsei", cex = 1.5,
pch = 16, ylim = c(4, 10))
# 1-st order
A <- cbind(rep(1, 5), x)
B <- y
cf <- lsei(A, B)$X
abline(coef = cf)
# 2-nd order
A <- cbind(A, x^2)
cf <- lsei(A, B)$X
curve(cf[1] + cf[2]*x + cf[3]*x^2, add = TRUE, lty = 2)
# 3-rd order
A <- cbind(A, x^3)
cf <- lsei(A, B)$X
curve(cf[1] + cf[2]*x + cf[3]*x^2 + cf[4]*x^3, add = TRUE, lty = 3)
# 4-th order
A <- cbind(A, x^4)
cf <- lsei(A, B)$X
curve(cf[1] + cf[2]*x + cf[3]*x^2 + cf[4]*x^3 + cf[5]*x^4,
add = TRUE, lty = 4)
legend("bottomleft", c("1st-order", "2nd-order","3rd-order","4th-order"),
lty = 1:4)
# ------------------------------------------------------------------------------
# example 2: equalities, approximate equalities and inequalities
# ------------------------------------------------------------------------------
A <- matrix(nrow = 4, ncol = 3,
data = c(3, 1, 2, 0, 2, 0, 0, 1, 1, 0, 2, 0))
B <- c(2, 1, 8, 3)
E <- c(0, 1, 0)
F <- 3
G <- matrix(nrow = 2, ncol = 3, byrow = TRUE,
data = c(-1, 2, 0, 1, 0, -1))
H <- c(-3, 2)
lsei(E = E, F = F, A = A, B = B, G = G, H = H)