# lsei: Least Squares with Equalities and Inequalities In limSolve: Solving Linear Inverse Models

## Description

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

 ```1 2 3``` ```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 <[email protected]>

## 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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46``` ```# ------------------------------------------------------------------------------ # 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) ```

### Example output

```Warning message:
In lsei(A, B) : No equalities - setting type = 2
Warning message:
In lsei(A, B) : No equalities - setting type = 2
Warning message:
In lsei(A, B) : No equalities - setting type = 2
Warning message:
In lsei(A, B) : No equalities - setting type = 2
\$X
[1]  1.2424242  3.0000000 -0.7575758

\$residualNorm
[1] 3.996803e-15

\$solutionNorm
[1] 98.06061

\$IsError
[1] FALSE

\$type
[1] "lsei"
```

limSolve documentation built on Aug. 14, 2017, 3:01 p.m.