LIN_SOL_LSQ
Solves a rectangular system of linear equations Ax ≅ b, in a least-squares sense. Using optional arguments, any of several related computations can be performed. These extra tasks include computing and saving the factorization of A using column and row pivoting, representing the determinant of A, computing the generalized inverse matrix A†, or computing the least-squares solution of
Ax ≅ b
or
ATy ≅ b,
given the factorization of A. An optional argument is provided for computing the following unscaled covariance matrix
Least-squares solutions, where the unknowns are non-negative or have simple bounds, can be computed with PARALLEL_NONNEGATIVE_LSQ and PARALLEL_BOUNDED_LSQ. These codes can be restricted to execute without MPI.
Required Arguments
A — Array of size m × n containing the matrix. (Input [/Output])
If the packaged option lin_sol_lsq_save_QR is used then the factorization of A is saved in A. For efficiency, the diagonal reciprocals of the matrix R are saved in the diagonal entries of A.
B — Array of size m × nb containing the right-hand side matrix. When using the option to solve adjoint systems ATx ≅ b, the size of b is n × nb. (Input [/Output])
If the packaged option lin_sol_lsq_save_QR is used then input B is used as work storage and is not saved.
X — Array of size n × nb containing the right-hand side matrix. When using the option to solve adjoint systems ATx ≅ b, the size of x is m × nb. (Output)
Optional Arguments
MROWS = m (Input)
Uses array A(1:m, 1:n) for the input matrix.
Default: m = size(A, 1)
NCOLS = n (Input)
Uses array A(1:m, 1:n) for the input matrix.
Default: n = size(A, 2)
NRHS = nb (Input)
Uses the array b(1:, 1:nb) for the input right-hand side matrix.
Default: nb = size(b, 2)
Note that b must be a rank-2 array.
pivots = pivots(:) (Output [/Input])
Integer array of size 2 * min(m, n) + 1 that contains the individual row followed by the column interchanges. The last array entry contains the approximate rank of A.
trans = trans(:) (Output [/Input])
Array of size 2 * min(m, n) that contains data for the construction of the orthogonal decomposition.
det = det(1:2) (Output)
Array of size 2 of the same type and kind as A for representing the products of the determinants of the matrices Q, P, and R. The determinant is represented by two numbers. The first is the base with the sign or complex angle of the result. The second is the exponent. When det(2) is within exponent range, the value of this expression is given by abs (det(1))**det(2) * (det(1))/abs(det(1)). If the matrix is not singular, abs(det(1)) = radix(det); otherwise, det(1) = 0, and det(2) = ‑ huge(abs(det(1))).
ainv = ainv(:,:) (Output)
Array with size n × m of the same type and kind as A(1:m, 1:n). It contains the generalized inverse matrix, A†.
cov = cov(:,:) (Output)
Array with size n × n of the same type and kind as A(1:m, 1:n). It contains the unscaled covariance matrix, C = (ATA)-1.
iopt = iopt(:) (Input)
Derived type array with the same precision as the input matrix; used for passing optional data to the routine. The options are as follows:
|
Packaged Options for lin_sol_lsq |
||
|
Option Prefix = ? |
Option Name |
Option Value |
|
s_, d_, c_, z_ |
lin_sol_lsq_set_small |
1 |
|
s_, d_, c_, z_ |
lin_sol_lsq_save_QR |
2 |
|
s_, d_, c_, z_ |
lin_sol_lsq_solve_A |
3 |
|
s_, d_, c_, z_ |
lin_sol_lsq_solve_ADJ |
4 |
|
s_, d_, c_, z_ |
lin_sol_lsq_no_row_pivoting |
5 |
|
s_, d_, c_, z_ |
lin_sol_lsq_no_col_pivoting |
6 |
|
s_, d_, c_, z_ |
lin_sol_lsq_scan_for_NaN |
7 |
|
s_, d_, c_, z_ |
lin_sol_lsq_no_sing_mess |
8 |
iopt(IO) = ?_options(?_lin_sol_lsq_set_small, Small)
Replaces with Small if a diagonal term of the matrix R is smaller in magnitude than the value Small. A solution is approximated based on this replacement in either case.
Default: the smallest number that can be reciprocated safely
iopt(IO) = ?_options(?_lin_sol_lsq_save_QR, ?_dummy)
Saves the factorization of A. Requires the optional arguments “pivots=” and “trans=” if the routine is used for solving further systems with the same matrix. This is the only case where the input arrays A and b are not saved. For efficiency, the diagonal reciprocals of the matrix R are saved in the diagonal entries of A.
iopt(IO) = ?_options(?_lin_sol_lsq_solve_A, ?_dummy)
Uses the factorization of A computed and saved to solve Ax = b.
iopt(IO) = ?_options(?_lin_sol_lsq_solve_ADJ, ?_dummy)
Uses the factorization of A computed and saved to solve ATx = b.
iopt(IO) = ?_options(?_lin_sol_lsq_no_row_pivoting, ?_dummy)
Does no row pivoting. The array pivots(:), if present, satisfies pivots(i) = i for i = 1, …, min (m, n).
iopt(IO) = ?_options(?_lin_sol_lsq_no_col_pivoting, ?_dummy)
Does no column pivoting. The array pivots(:), if present, satisfies pivots(i + min (m, n)) = i for i = 1, …, min (m, n).
iopt(IO) = ?_options(?_lin_sol_lsq_scan_for_NaN, ?_dummy)
Examines each input array entry to find the first value such that
isNaN(a(i,j)) .or. isNan(b(i,j)) == .true.
See the isNaN() function, Chapter 10.
Default: Does not scan for NaNs
iopt(IO) = ?_options(?_lin_sol_lsq_no_sing_mess, ?_dummy)
Do not print an error message when A is singular or k < min(m, n).
FORTRAN 90 Interface
Generic: CALL LIN_SOL_LSQ (A, B, X [, …])
Specific: The specific interface names are S_LIN_SOL_LSQ, D_LIN_SOL_LSQ, C_LIN_SOL_LSQ, and Z_LIN_SOL_LSQ.
Description
Routine LIN_SOL_LSQ solves a rectangular system of linear algebraic equations in a least-squares sense. It computes the decomposition of A using an orthogonal factorization. This decomposition has the form
where the matrices Q and P are products of elementary orthogonal and permutation matrices. The matrix R is k × k, where k is the approximate rank of A. This value is determined by the value of the parameter Small. See Golub and Van Loan (1989, Chapter 5.4) for further details. Note that the use of both row and column pivoting is nonstandard, but the routine defaults to this choice for enhanced reliability.
Fatal and Terminal Error Messages
See the messages.gls file for error messages for LIN_SOL_LSQ. These error messages are numbered 241–256; 261–276; 281–296; 301–316.
Examples
Example 1: Solving a Linear Least-squares System
This example solves a linear least-squares system Cx ≅ d, where
is a real matrix with m > n. The least-squares problem is derived from polynomial data fitting to the function
using a discrete set of values in the interval –1 ≤ x ≤ 1. The polynomial is represented as the series
where the Ti(x) are Chebyshev polynomials. It is natural for the problem matrix and solution to have a column or entry corresponding to the subscript zero, which is used in this code. Also, see operator_ex09, supplied with the product examples.
use lin_sol_lsq_int
use rand_gen_int
use error_option_packet
implicit none
! This is Example 1 for LIN_SOL_LSQ.
integer i
integer, parameter :: m=128, n=8
real(kind(1d0)), parameter :: one=1d0, zero=0d0
real(kind(1d0)) A(m,0:n), c(0:n,1), pi_over_2, x(m), y(m,1), &
u(m), v(m), w(m), delta_x
! Generate a random grid of points.
call rand_gen(x)
! Transform points to the interval -1,1.
x = x*2 - one
! Compute the constant 'PI/2'.
pi_over_2 = atan(one)*2
! Generate known function data on the grid.
y(1:m,1) = exp(x) + cos(pi_over_2*x)
! Fill in the least-squares matrix for the Chebyshev polynomials.
A(:,0) = one; A(:,1) = x
do i=2, n
A(:,i) = 2*x*A(:,i-1) - A(:,i-2)
end do
! Solve for the series coefficients.
call lin_sol_lsq(A, y, c)
! Generate an equally spaced grid on the interval.
delta_x = 2/real(m-1,kind(one))
do i=1, m
x(i) = -one + (i-1)*delta_x
end do
! Evaluate residuals using backward recurrence formulas.
u = zero
v = zero
do i=n, 0, -1
w = 2*x*u - v + c(i,1)
v = u
u = w
end do
y(1:m,1) = exp(x) + cos(pi_over_2*x) - (u-x*v
! Check that n+1 sign changes in the residual curve occur.
x = one
x = sign(x,y(1:m,1))
if (count(x(1:m-1) /= x(2:m)) >= n+1) then
write (*,*) 'Example 1 for LIN_SOL_LSQ is correct.'
end if
end
Example 1 for LIN_SOL_LSQ is correct.
Example 2: System Solving with the Generalized Inverse
This example solves the same form of the system as Example 1. In this case, the grid of evaluation points is equally spaced. The coefficients are computed using the “smoothing formulas” by rows of the generalized inverse matrix, A†, computed using the optional argument “ainv=”. Thus, the coefficients are given by the matrix-vector product c = (A†) y, where y is the vector of values of the function y(x) evaluated at the grid of points. Also, see operator_ex10, supplied with the product examples.
use lin_sol_lsq_int
implicit none
! This is Example 2 for LIN_SOL_LSQ.
integer i
integer, parameter :: m=128, n=8
real(kind(1d0)), parameter :: one=1.0d0, zero=0.0d0
real(kind(1d0)) a(m,0:n), c(0:n,1), pi_over_2, x(m), y(m,1), &
u(m), v(m), w(m), delta_x, inv(0:n, m)
! Generate an array of equally spaced points on the interval -1,1.
delta_x = 2/real(m-1,kind(one))
do i=1, m
x(i) = -one + (i-1)*delta_x
end do
! Compute the constant 'PI/2'.
pi_over_2 = atan(one)*2
! Compute data values on the grid.
y(1:m,1) = exp(x) + cos(pi_over_2*x)
! Fill in the least-squares matrix for the Chebyshev polynomials.
a(:,0) = one
a(:,1) = x
do i=2, n
a(:,i) = 2*x*a(:,i-1) - a(:,i-2)
end do
! Compute the generalized inverse of the least-squares matrix.
call lin_sol_lsq(a, y, c, nrhs=0, ainv=inv)
! Compute the series coefficients using the generalized inverse
! as 'smoothing formulas.'
c(0:n,1) = matmul(inv(0:n,1:m),y(1:m,1))
! Evaluate residuals using backward recurrence formulas.
u = zero
v = zero
do i=n, 0, -1
w = 2*x*u - v + c(i,1)
v = u
u = w
end do
y(1:m,1) = exp(x) + cos(pi_over_2*x) - (u-x*v)
! Check that n+2 sign changes in the residual curve occur.
! (This test will fail when n is larger.)
x = one
x = sign(x,y(1:m,1))
if (count(x(1:m-1) /= x(2:m)) == n+2) then
write (*,*) 'Example 2 for LIN_SOL_LSQ is correct.'
end if
end
Output
Example 2 for LIN_SOL_LSQ is correct.
Example 3: Two-Dimensional Data Fitting
This example illustrates the use of radial-basis functions to least-squares fit arbitrarily spaced data points. Let m data values {yi} be given at points in the unit square, {pi}. Each pi is a pair of real values. Then, n points {qj} are chosen on the unit square. A series of radial-basis functions is used to represent the data,
where δ2 is a parameter. This example uses δ2 = 1, but either larger or smaller values can give a better approximation for user problems. The coefficients {cj} are obtained by solving the following m × n linear least-squares problem:
This example illustrates an effective use of Fortran 90 array operations to eliminate many details required to build the matrix and right-hand side for the {cj} . For this example, the two sets of points {pi} and {qj} are chosen randomly. The values {yj} are computed from the following formula:
The residual function
is computed at an N × N square grid of equally spaced points on the unit square. The magnitude of r(p) may be larger at certain points on this grid than the residuals at the given points, {pi}. Also, see operator_ex11, supplied with the product examples.
use lin_sol_lsq_int
use rand_gen_int
implicit none
! This is Example 3 for LIN_SOL_LSQ.
integer i, j
integer, parameter :: m=128, n=32, k=2, n_eval=16
real(kind(1d0)), parameter :: one=1.0d0, delta_sqr=1.0d0
real(kind(1d0)) a(m,n), b(m,1), c(n,1), p(k,m), q(k,n), &
x(k*m), y(k*n), t(k,m,n), res(n_eval,n_eval), &
w(n_eval), delta
! Generate a random set of data points in k=2 space.
call rand_gen(x)
p = reshape(x,(/k,m/))
! Generate a random set of center points in k-space.
call rand_gen(y)
q = reshape(y,(/k,n/))
! Compute the coefficient matrix for the least-squares system.
t = spread(p,3,n)
do j=1, n
t(1:,:,j) = t(1:,:,j) - spread(q(1:,j),2,m)
end do
a = sqrt(sum(t**2,dim=1) + delta_sqr)
! Compute the right hand side of data values.
b(1:,1) = exp(-sum(p**2,dim=1))
! Compute the solution.
call lin_sol_lsq(a, b, c)
! Check the results.
if (sum(abs(matmul(transpose(a),b-matmul(a,c))))/sum(abs(a)) &
<= sqrt(epsilon(one))) then
write (*,*) 'Example 3 for LIN_SOL_LSQ is correct.'
end if
! Evaluate residuals, known function - approximation at a square
! grid of points. (This evaluation is only for k=2.)
delta = one/real(n_eval-1,kind(one))
do i=1, n_eval
w(i) = (i-1)*delta
end do
res = exp(-(spread(w,1,n_eval)**2 + spread(w,2,n_eval)**2))
do j=1, n
res = res - c(j,1)*sqrt((spread(w,1,n_eval) - q(1,j))**2 + &
(spread(w,2,n_eval) - q(2,j))**2 + delta_sqr)
end do
end
Output
Example 3 for LIN_SOL_LSQ is correct.
Example 4: Least-squares with an Equality Constraint
This example solves a least-squares system Ax ≅ b with the constraint that the solution values have a sum equal to the value 1. To solve this system, one heavily weighted row vector and right-hand side component is added to the system corresponding to this constraint. Note that the weight used is
where Ɛ is the machine precision, but any larger value can be used. The fact that lin_sol_lsq performs row pivoting in this case is critical for obtaining an accurate solution to the constrained problem solved using weighting. See Golub and Van Loan (1989, Chapter 12) for more information about this method. Also, see operator_ex12, supplied with the product examples.
use lin_sol_lsq_int
use rand_gen_int
implicit none
! This is Example 4 for LIN_SOL_LSQ.
integer, parameter :: m=64, n=32
real(kind(1e0)), parameter :: one=1.0e0
real(kind(1e0)) :: a(m+1,n), b(m+1,1), x(n,1), y(m*n)
! Generate a random matrix.
call rand_gen(y)
a(1:m,1:n) = reshape(y,(/m,n
! Generate a random right hand side.
call rand_gen(b(1:m,1))
! Heavily weight desired constraint. All variables sum to one.
a(m+1,1:n) = one/sqrt(epsilon(one))
b(m+1,1) = one/sqrt(epsilon(one))
call lin_sol_lsq(a, b, x)
if (abs(sum(x) - one)/sum(abs(x)) <= &
sqrt(epsilon(one))) then
write (*,*) 'Example 4 for LIN_SOL_LSQ is correct.'
end if
end
Output
Example 4 for LIN_SOL_LSQ is correct.
