"Hybrid" LS-LP model pseudoinverse-based (SVD-based) solving algorithm
Project description
The algorithm solves “hybrid” least squares linear programming (LS-LP) problems with the help of the Moore-Penrose inverse (pseudoinverse), calculated using singular value decomposition (SVD), with emphasis on estimation of non-typical constrained OLS (cOLS), Relationship Matrix (RM) [1], and custom (user-defined) cases. The pseudoinverse offers a unique solution and may be the best linear unbiased estimator (BLUE) for a group of problems under certain conditions [2].
Mechanism:
The problem is written as a matrix equation a @ x = b where a consists of coefficients for CONSTRAINTS and for SLACK VARIABLES (the upper part) as well as for MODEL (the lower part) as illustrated in Figure 1. Each part of a can be omitted to accommodate a special case:
cOLS problems require no case-specific CONSTRAINTS;
RM problems require case-specific CONSTRAINTS, no problem CONSTRAINTS, and an optional MODEL;
SLACK VARIABLES are non-zero only for inequality constraints and are omitted if problems don’t include any;
…
a |
b |
---|---|
CONSTRAINTS (PROBLEM + CASE-SPECIFIC) | SLACK VARIABLES |
CONSTRAINTS |
MODEL |
MODEL |
Source: self-prepared
Parameters:
# LS-LP type and input
lp : str or unicode, optional. ‘cOLS’ (Constrained OLS), ‘RM’ (Relationship Matrix), or ‘custom’ (user-defined case).
lhs : sequence of array_like. Left-hand side of the problem (NB 3D iterable for lp=‘cOLS’).
rhs : sequence of array_like. Right-hand side of the problem (NB row sums first for lp=‘RM’).
svar : array_like. Slack variables.
# a and estimation
diag : bool, optional. Diagonal of M x N (NB for lp=‘RM’).
rcond : float, optional. Cut-off ratio for small singular values of a. For the purposes of rank determination, singular values are treated as zero if they are smaller than rcond times the largest singular value of a in numpy.linalg.lstsq(a, b, rcond=rcond).
Returns:
x : {(N,), (N, K)} ndarray. Least-squares solution.
mse : float. Sums of squared residuals: Squared Euclidean 2-norm for each column in `b - a @ x. If the rank of a is < N or M <= N, a problem-specific formula is applied. NB Empty for lp=‘custom’ (must be calculated by user).
a : (ndarray, int) tuple. Matrix a, Rank of a.
s : (min(M, N),) ndarray. Singular values of a.
Raises:
LPpinvError If the LS-LP problem is incorrectly specified.
LinAlgError If computation does not converge.
Examples:
import numpy as np import lppinv as lp # cOLS problem lp.solve( lp='cOLS', lhs=([[0, 1, 2, 3], [0, 5, 2, 8], np.ones(4)],[[0, 1, 2, 3], [0, 5, 2, 8], np.ones(4)]), rhs=[np.zeros(4), [-1, 0.2, 0.9, 2.1]], svar=[-1, 1, -1, 0] ) # RM problem lp.solve( lp='RM', rhs=[[4, 5, 3, 4, 6], [1, 2], [0, 0]], diag=True ) # custom problem lp.solve( lp='custom', lhs=np.vstack([[0, 1, 1], [1, 0, 1], [1, 1, 0]]), rhs=np.array([2, 3, 9]).T, svar=[-1, 0, 0] )
Notes:
[1] Relationship Matrix (RM) of size M x N is a formal model of relationships between M and N elements in a system. For example,
an input-output table (IOT) is a type of RM where M = N and the elements are industries;
a matrix of trade / investment / etc. is a type of RM where M = N and the elements are countries or (macro)regions in which diagonal elements can, in some cases, all be equal to zero.
[2] For example, consult Albert, A., 1972. Regression And The Moore-Penrose Pseudoinverse. New York: Academic Press. Chapter VII.
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.