lin_least_squares_gen
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 imsl_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)
float tol (Input)
Nonnegative tolerance used to determine the subset of columns of A to be included in the solution.
Default: tol = sqrt (imsl_amach(4))
int *kbasis (Input/Output)
Integer containing the number of columns used in the solution.
kbasis = k if
∣rk+1,k+1∣ < |
tol∣*∣r1,1∣. 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)
float **p_qraux (Input/Output)
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.
float **p_qr (Input/Output)
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)
float qraux[] (Input/Output)
A user-allocated array of size n containing the scalars τk of the Householder transformations in the first min (m, n) positions.
float qr[] (Input/Output)
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_USER, 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};
int 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};
int 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>
#include <stdio.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};
int 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 imsl_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 #. |
