linLeastSquaresGen¶
Solves a linear least-squares problem \(Ax = b\). Using optional arguments, the QR factorization of A, \(AP = QR\), and the solve step based on this factorization can be computed.
Synopsis¶
linLeastSquaresGen (a, b)
Required Arguments¶
- float
a[[]](Input) - Array of size m × n containing the matrix.
- float
b[](Input) - Array of size m containing the right-hand side.
Return Value¶
If no optional arguments are used, function linLeastSquaresGen returns
the solution x of the linear least-squares problem \(Ax = b\). If no value
can be computed, then None is returned.
Optional Arguments¶
basis, float tol, int kbasis (Input, Input/Output)
tol(Input)Nonnegative tolerance used to determine the subset of columns of A to be included in the solution.
Default: tol = sqrt (
machine(4))kbasis(Input/Output)Integer containing the number of columns used in the solution.
kbasis= k if \(|r_{k+1,k+1}| < |\textit{tol}| * |r_{1,1}|\). For more information on the use of this option, see the Description section.Default:
kbasis= min (m, n)
residual(Output)- An array of size m containing the residual vector b − Ax.
factor, pQraux, pQr (Output)
- float
pQraux(Input/Output)- An array of size n containing the scalars \(τ_k\) of the Householder transformations in the first min (m, n) positions.
- float
pQr(Input/Output)- An array of size m × n containing the Householder transformations that define the decomposition. The strictly lower-triangular part of this array contains the information to construct Q, and the upper-triangular part contains R.
These parameters are “Input” if
solveis specified; “Output” otherwise.
q(Output)- An array of size m × m containing the orthogonal matrix of the factorization.
pivot, intpivot[](Input/Output)Array of size n containing the desired variable order and usage information. The argument is used with
factorOnlyorsolveOnly.On input, if
pivot[k − 1] > 0, then column k of A is an initial column. Ifpivot[k − 1] = 0, then the column of A is a free column and can be interchanged in the column pivoting. Ifpivot[k − 1] < 0, then column k of A is a final column. If all columns are specified as initial (or final) columns, then no pivoting is performed. (The permutation matrix P is the identity matrix in this case.)On output,
pivot[k − 1] contains the index of the column of the original matrix that has been interchanged into column k.Default:
pivot[k − 1] = 0, k = 1, …, nfactorOnly- Compute just the QR factorization of the matrix AP with the
permutation matrix P defined by
pivotand by further pivoting involving free columns. IffactorOnlyis used, the additional argumentspivotandfactorare required. In that case, the required argumentbis ignored, and the returned value of the function isNone. solveOnly- Compute the solution to the least-squares problem \(Ax = b\) given the QR factorization previously computed by this function.
Description¶
The function linLeastSquaresGen solves a system of linear least-squares
problems \(Ax = b\) with column pivoting. It computes a QR factorization
of the matrix AP, where P is the permutation matrix defined by the
pivoting, and computes the smallest integer k satisfying
\(|r_{k+1,k+1}| < ∣\mathit{tol}∣ * |r_{1,1}|\) to the output variable
kbasis. Householder transformations
\(k = 1, …, min (m−1, n)\) are used to compute the factorization. The
decomposition is computed in the form \(Q_{min(m-1,n)} \ldots Q_1AP =
R\), so \(AP = QR\) where \(Q = Q_1 \ldots Q_{min(m-1,n)}\). Since each
Householder vector \(u_k\) has zeros in the first k − 1 entries, it is
stored as part of column k of qr. The upper-trapezoidal matrix R is
stored in the upper-trapezoidal part of the first min (m, n) rows of
qr. The solution x to the least-squares problem is computed by solving
the upper-triangular system of linear equations \(R(1:k,1:k)y(1:k) =
(Q^Tb)(1:k)\) with k = kbasis. The solution is completed by setting
\(y(k+1 : n)\) to zero and rearranging the variables, \(x = Py\).
When factorOnly is specified, the function computes the QR
factorization of AP with P defined by the input pivot and by column
pivoting among ‘‘free’’ columns. Before the factorization, initial columns
are moved to the beginning of the array a and the final columns to the
end. Both initial and final columns are not permuted further during the
computation. Just the free columns are moved.
If solveOnly is specified, then the function computes the least-squares
solution to \(Ax = b\) given the QR factorization previously defined.
There are kbasis columns used in the solution. Hence, in the case that
all columns are free, x is computed as described in the default case.
Examples¶
Example 1¶
This example illustrates the least-squares solution of four linear equations
in three unknowns using column pivoting. The problem is equivalent to
least-squares quadratic polynomial fitting to four data values. Write the
polynomial as \(p(t) = x_1 + tx_2 + t^2 x_3\) and the data pairs
(\(t_i\), \(b_i\)), \(t_i = 2i, i = 1, 2, 3, 4\). The solution to
\(Ax = b\) is returned by the function linLeastSquaresGen.
from numpy import *
from pyimsl.math.linLeastSquaresGen import linLeastSquaresGen
from pyimsl.math.writeMatrix import writeMatrix
a = array([[1.0, 2.0, 4.0],
[1.0, 4.0, 16.0],
[1.0, 6.0, 36.0],
[1.0, 8.0, 64.0]])
b = array([4.999, 9.001, 12.999, 17.001])
# Solve Ax=b for x
x = linLeastSquaresGen(a, b)
writeMatrix("Solution vector", x)
Output¶
Solution vector
1 2 3
0.999 2.000 -0.000
Example 2¶
This example uses the same coefficient matrix A as in the initial example.
It computes the QR factorization of A with column pivoting. The final
and free columns are specified by pivot and the column pivoting is done
only among the free columns.
from numpy import *
from pyimsl.math.linLeastSquaresGen import linLeastSquaresGen
from pyimsl.math.writeMatrix import writeMatrix
a = array([[1.0, 2.0, 4.0],
[1.0, 4.0, 16.0],
[1.0, 6.0, 36.0],
[1.0, 8.0, 64.0]])
b = []
pvt = [0, 0, -1]
factor = {}
p_q = []
# Compute the QR factorization of A with partial column pivoting
x = linLeastSquaresGen(a, b,
pivot=pvt,
factor=factor,
q=p_q,
factorOnly=True)
writeMatrix("The matrix Q", p_q)
writeMatrix("The matrix R", factor['pQr'], printUpper=True)
writeMatrix("The Pivot Sequence ", pvt)
Output¶
The matrix Q
1 2 3 4
1 -0.1826 -0.8165 0.5000 -0.2236
2 -0.3651 -0.4082 -0.5000 0.6708
3 -0.5477 0.0000 -0.5000 -0.6708
4 -0.7303 0.4082 0.5000 0.2236
The matrix R
1 2 3
1 -10.95 -1.83 -73.03
2 -0.82 16.33
3 8.00
The Pivot Sequence
1 2 3
2 1 3
Example 3¶
This example computes the QR factorization with column pivoting for the matrix A of the initial example. It computes the least-squares solutions to \(Ax = b_i\) for \(i = 1, 2, 3\).
from __future__ import print_function
from numpy import *
from pyimsl.math.linLeastSquaresGen import linLeastSquaresGen
from pyimsl.math.writeMatrix import writeMatrix
a = array([[1.0, 2.0, 4.0],
[1.0, 4.0, 16.0],
[1.0, 6.0, 36.0],
[1.0, 8.0, 64.0]])
b = array([4.999, 9.001, 12.999, 17.001,
2.0, 3.142, 5.11, 0.0,
1.34, 8.112, 3.76, 10.99])
m = 4
n = 3
k = 3
pvt = [0, 0, 0]
p_qraux = []
p_qr = []
basis = {'tol': 1.e-4}
factor = {}
# Factor A with the given pvt setting all variables to be free
linLeastSquaresGen(a, b[0:4],
basis=basis,
pivot=pvt,
factor=factor,
factorOnly=True)
# Print some factorization information
print("Number of Columns in the base %2d" % (basis['kbasis']))
writeMatrix("Upper triangular R Matrix", factor['pQr'],
printUpper=True)
writeMatrix("The output column order ", pvt)
# Solve Ax = b for each x given the factorization
p_res = []
for i in range(0, k):
x = linLeastSquaresGen(a, b[(i * m):(m * (i + 1))],
basis=basis,
pivot=pvt,
factor=factor,
residual=p_res,
solveOnly=True)
# Print right-hand side, b and solution, x
writeMatrix("Right side,B", b[(i * m):(m * (i + 1))])
writeMatrix("Solution, x ", x)
# Print residuals, b - Ax
writeMatrix("Residual, b - Ax ", p_res)
Output¶
Number of Columns in the base 3
Upper triangular R Matrix
1 2 3
1 -75.26 -10.63 -1.59
2 -2.65 -1.15
3 0.36
The output column order
1 2 3
3 2 1
Right side,B
1 2 3 4
5 9 13 17
Solution, x
1 2 3
0.999 2.000 -0.000
Residual, b - Ax
1 2 3 4
-0.0004 0.0012 -0.0012 0.0004
Right side,B
1 2 3 4
2.000 3.142 5.110 0.000
Solution, x
1 2 3
-4.244 3.706 -0.391
Residual, b - Ax
1 2 3 4
0.395 -1.186 1.186 -0.395
Right side,B
1 2 3 4
1.34 8.11 3.76 10.99
Solution, x
1 2 3
0.4735 0.9436 0.0286
Residual, b - Ax
1 2 3 4
-1.135 3.406 -3.406 1.135
Fatal Errors¶
IMSL_SINGULAR_TRI_MATRIX |
The input triangular matrix is singular. The index of the first zero diagonal term is #. |