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
#include <imsl.h>
float *imsl_f_lin_least_squares_gen (int m, int n, float a[], float b[], …, 0)
The type double procedure is imsl_d_lin_least_squares_gen.
Required Arguments
int m
(Input)
Number of rows in the matrix.
int n
(Input)
Number of columns in the matrix.
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 imsl_f_lin_least_squares_gen returns a pointer to the solution x of the linear least-squares problem Ax = b. To release this space, use free. If no value can be computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
float
*imsl_f_lin_least_squares_gen (int
m, int n, float
a[],
float b[],
IMSL_A_COL_DIM,
int
a_col_dim,
IMSL_RETURN_USER, float
x[],
IMSL_BASIS, float
tol,
int
*kbasis,
IMSL_RESIDUAL, float
**p_res,
IMSL_RESIDUAL_USER, float
res[],
IMSL_FACTOR, float
**p_qraux,
float
**p_qr,
IMSL_FACTOR_USER, float
qraux[],
float
qr[],
IMSL_FAC_COL_DIM, int
qr_col_dim,
IMSL_Q, float
**p_q,
IMSL_Q_USER, float
q[],
IMSL_Q_COL_DIM, int
q_col_dim,
IMSL_PIVOT, int
pvt[],
IMSL_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
0)
Optional Arguments
IMSL_A_COL_DIM, int a_col_dim
(Input)
The column dimension of the array a.
Default: a_col_dim = n
IMSL_RETURN_USER, float x[]
(Output)
A user-allocated array of size n containing the least-squares
solution x. If IMSL_RETURN_USER is
used, the return value of the function is a pointer to the array x.
IMSL_BASIS, float tol, int *kbasis (Input, Input/Output)
tol: Nonnegative tolerance used to determine the subset of columns of A to be included in the solution.
Default: tol = sqrt (imsl_amach(4))
kbasis: Integer containing the number of columns used in the solution.
kbasis = k if |rk+1,+1| < |tol|*|r1.1| and |ri.i|³ tol*|r1.1| for i = 1, 2, …, k. For more information on the use of this option, see “Description” section.
Default: kbasis = min (m, n)
IMSL_RESIDUAL, float **p_res
(Output)
The address of a pointer to an array of size m containing the
residual vector b − Ax. On return,
the necessary space is allocated by the function. Typically, float
*p_res is
declared, and &p_res is used as
an argument.
IMSL_RESIDUAL_USER, float res[]
(Output)
A user-allocated array of size m containing the residual
vector b − Ax.
IMSL_FACTOR, float **p_qraux, float **p_qr (Output)
**p_qraux: The address of a pointer qraux to an array of size n containing the scalars τk of the Householder transformations in the first min (m, n) positions. On return, the necessary space is allocated by the function. Typically, float *qraux is declared, and &qraux is used as an argument.
**p_qr: The address of a pointer to 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. On return, the necessary space is allocated by the function. Typically, float *qr is declared, and &qr is used as an argument.
IMSL_FACTOR_USER, float qraux[], float qr[] (Input /Output)
qraux[]: A user-allocated array of size n containing the scalars τk of the Householder transformations in the first min (m, n) positions.
qr[]: A user-allocated 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. The upper-triangular part contains R. If the data in a is not needed, qr can share the same storage locations as a by using a instead of the separate argument qr.
These parameters are “Input” if IMSL_SOLVE is specified; “Output” otherwise.
IMSL_FAC_COL_DIM, int
qr_col_dim (Input)
The column dimension of the array
containing QR factorization.
Default: qr_col_dim = n
IMSL_Q, float **p_q
(Output)
The address of a pointer to an array of size
m × m containing the orthogonal matrix of the
factorization. On return, the necessary space is allocated by the function.
Typically, float
*q is declared, and &q is used as an
argument.
IMSL_Q_USER, float q[]
(Output)
A user-allocated array of size m × m
containing the orthogonal matrix Q of the QR factorization.
IMSL_Q_COL_DIM, int q_col_dim
(Input)
The column dimension of the array containing the Q matrix of
the factorization.
Default: q_col_dim = m
IMSL_PIVOT, int pvt[]
(Input/Output)
Array of size n containing the desired variable order
and usage information. The argument is used with IMSL_FACTOR_ONLY or
IMSL_SOLVE_ONLY.
On input, if pvt [k − 1] > 0, then column k of A is an initial column. If pvt [k − 1] = 0, then the column of A is a free column and can be interchanged in the column pivoting. If pvt [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, pvt [k − 1] contains the index of the column of the original matrix that has been interchanged into column k.
Default: pvt [k − 1] = 0, k = 1, …, n
IMSL_FACTOR_ONLY
Compute
just the QR factorization of the matrix AP with the permutation
matrix P defined by pvt and by further
pivoting involving free columns. If IMSL_FACTOR_ONLY is
used, the additional arguments IMSL_PIVOT and IMSL_FACTOR are
required. In that case, the required argument b is ignored, and the
returned value of the function is NULL.
IMSL_SOLVE_ONLY
Compute
the solution to the least-squares problem Ax = b given
the QR factorization previously computed by this function. If IMSL_SOLVE_ONLY is
used, arguments IMSL_FACTOR, IMSL_PIVOT, and IMSL_BASIS are
required, and the required argument a is ignored.
Description
The function imsl_f_lin_least_squares_gen 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 |rk+1,k+1| < |tol|*|r1,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 Qmin(m-1, n)…Q1AP = R,
so AP = QR where Q = Q1…Qmin(m-1, n). Since
each Householder vector uk 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) = (QTb) (1:k)
with k = kbasis.
The solution is completed by setting
y(k + 1 : n) to
zero and rearranging the variables, x = Py.
When IMSL_FACTOR_ONLY is specified, the function computes the QR factorization of AP with P defined by the input pvt 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 IMSL_SOLVE_ONLY 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) = x1 + tx2 + t2x3 and the data pairs (ti, bi), ti = 2i, i = 1, 2, 3, 4. A pointer to the solution to Ax = b is returned by the function imsl_f_lin_least_squares_gen.
#include <imsl.h>
float
a[] = {1.0, 2.0,
4.0,
1.0, 4.0,
16.0,
1.0, 6.0,
36.0,
1.0, 8.0, 64.0};
float b[] = {4.999, 9.001,
12.999, 17.001};
main()
{
int m = 4, n =
3;
float
*x;
/* Solve Ax = b for x */
x =
imsl_f_lin_least_squares_gen (m, n, a, b,
0);
/* Print x */
imsl_f_write_matrix ("Solution vector", 1,
n, x, 0);
}
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 pvt and the column pivoting is done only among the free columns.
#include
<imsl.h>
float a[] = {1.0, 2.0,
4.0,
1.0, 4.0,
16.0,
1.0, 6.0,
36.0,
1.0, 8.0, 64.0};
int pvt[] = {0, 0,
-1};
main()
{
int m = 4, n =
3;
float
*x, *b;
float *p_qraux,
*p_qr;
float
*p_q;
/* Compute the QR factorization
*/
/* of A with partial column
*/
/* pivoting */
x = imsl_f_lin_least_squares_gen (m, n, a,
b,
IMSL_PIVOT,
pvt,
IMSL_FACTOR, &p_qraux,
&p_qr,
IMSL_Q,
&p_q,
IMSL_FACTOR_ONLY,
0);
/* Print Q */
imsl_f_write_matrix ("The matrix Q", m, m,
p_q,
0);
/*
Print R */
imsl_f_write_matrix ("The matrix R", m, n,
p_qr,
IMSL_PRINT_UPPER,
0);
/* Print pivots */
imsl_i_write_matrix ("The Pivot
Sequence", 1, n, pvt, 0);
}
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 = bi for i = 1, 2, 3.
#include <imsl.h>
float
a[] = {1.0, 2.0,
4.0,
1.0, 4.0,
16.0,
1.0, 6.0,
36.0,
1.0, 8.0, 64.0};
float b[] = {4.999, 9.001,
12.999,
17.001,
2.0, 3.142, 5.11,
0.0,
1.34, 8.112, 3.76,
10.99};
int pvt[] = {0, 0,
0};
main()
{
int m = 4, n =
3;
int i, k =
3;
float
*p_qraux, *p_qr;
float tol =
1.e-4;
int
*kbasis;
float *x,
*p_res;
/* Factor A with the given pvt
*/
/* setting all variables to
*/
/* be free */
imsl_f_lin_least_squares_gen (m, n, a, b,
IMSL_BASIS,
tol, &kbasis,
IMSL_PIVOT,
pvt,
IMSL_FACTOR,
&p_qraux, &p_qr,
IMSL_FACTOR_ONLY,
0);
/* Print some factorization
*/
/* information*/
printf("Number of Columns in the
base\n%2d", kbasis);
imsl_f_write_matrix ("Upper
triangular R Matrix", m, n, p_qr,
IMSL_PRINT_UPPER,
0);
imsl_i_write_matrix ("The output column order ",
1, n, pvt,
0);
/* Solve Ax = b for each x
*/
/* given the factorization */
for ( i = 0; i < k; i++)
{
x =
imsl_f_lin_least_squares_gen (m, n, a, &b[i*m],
IMSL_BASIS, tol,
&kbasis,
IMSL_PIVOT,
pvt,
IMSL_FACTOR_USER, p_qraux,
p_qr,
IMSL_RESIDUAL, &p_res,
IMSL_SOLVE_ONLY,
0);
/* Print right-hand side, b
*/
/* and solution, x */
imsl_f_write_matrix ("Right-hand side, b ", 1, m,
&b[i*m], 0);
imsl_f_write_matrix ("Solution, x ", 1, n, x,
0);
/* Print residuals, b - Ax
*/
imsl_f_write_matrix
("Residual, b - Ax ", 1, m, p_res,
0);
}
}
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-hand 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-hand 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-hand side, b
1
2
3
4
1.34
8.11
3.76 10.99
Solution, x
1
2
3
0.4735
0.9437
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 #.
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |