Module_LLSQ_Procedures¶

This file is part of statpack.

statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)

MODULE EXPORTING SUBROUTINES AND FUNCTIONS FOR SOLVING LINEAR LEAST SQUARES PROBLEMS.

LATEST REVISION : 12/03/2024

`function solve_llsq ( a, b, krank, tol, min_norm )`¶

Purpose¶

SOLVE_LLSQ computes a solution to a real linear least squares problem:

Minimize 2-norm || B - A * X ||

using an orthogonal factorization with columns pivoting of A. A is a m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted. Here, B is a m-element right hand side vector and X is a n-element solution vector.

The function returns the n-element solution vector X.

A and B are not overwritten by SOLVE_LLSQ.

Arguments¶

A (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-n coefficient matrix A.

B (INPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

The shape of B must verify:

size( B ) = size( A, 1 ) = m .

KRANK (INPUT, OPTIONAL) integer(i4b)

On entry, KRANK=k implies that the first k columns of A are to be forced into the basis, pivoting is performed on the last n-k columns of A.

When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be full rank.

If KRANK is absent or is <=0, pivoting is done on all columns of A. This is appropriate if A may not be of full rank (i.e. certain columns of A are linear combinations of other columns).

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of A, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of A whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.

TOL must not be greater or equal to 1 or less or equal than 0, otherwise the numerical rank of A is determined and the calculations to determine the condition number are not performed. If TOL is absent, the numerical rank of A is determined.

MIN_NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

MIN_NORM=true, the minimun 2-norm solution is computed.

MIN_NORM=false or if MIN_NORM is absent, a solution is computed such that if the j-th column of A is omitted from the basis, X[j] is set to zero.

Further Details¶

The routine first computes a QR factorization with (partial) column pivoting on option (see below):

A * P = Q * R

, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.

R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, arank, is the effective rank of A.

This leads to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.

If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

A * P = Q * T * Z

, where Z is a n-by-n orthogonal matrix and T has the form:

[ T11 T12 ]

[ T21 T22 ]

, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a arank-by-arank upper triangular matrix.

The minimum 2-norm solution is then

X = [ P * Z’ ](:,:arank) * [ inv(T11) * Q1’ * B ]

where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first arank columns of Q.

If MIN_NORM is absent or has the value false, a solution is computed as

X = P(:,:arank) * [ inv(R11) * Q1’ * B ]

where inv(R11) is the inverse of R11 and Q1 consists of the first arank columns of Q. In this case, if the j-th column of A is omitted from the basis, X[j] is set to zero.

On input, if KRANK is present and KRANK=k, the first k columns of A are to be forced into the basis. Pivoting is performed on the last n-k columns of A.

When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be of full rank.

If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of A.

TOL is an optional argument such that 0<TOL<1. If TOL is not specified, or is outside ]0,1[, the calculations to determine the condition number of A are not performed and crude tests on R(j,j) are performed to determine the numerical rank of A. If TOL is present and is in ]0,1[, the calculations to determine the condition number are performed.

If it is possible that A may not be full rank (i.e., certain columns of A are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of A so that all elements are about the same order of magnitude.

On exit, if A or B are empty, the function returns a n-element vector filled with nan() value.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`function solve_llsq ( a, b, krank, tol, min_norm )`¶

Purpose¶

SOLVE_LLSQ computes solutions to real linear least squares problems of the form:

Minimize 2-norm|| B - A * X ||

using an orthogonal factorization with columns pivoting of A. A is a m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted.

Several right hand side vectors b can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B.

The function returns the n-by-nb solution matrix X.

A and B are not overwritten by SOLVE_LLSQ.

Arguments¶

A (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-n coefficient matrix A.

B (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B.

The shape of B must verify:

size( B, 1 ) = size( A, 1 ) = m .

KRANK (INPUT, OPTIONAL) integer(i4b)

On entry, KRANK=k implies that the first k columns of A are to be forced into the basis, pivoting is performed on the last n-k columns of A.

When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be full rank.

If KRANK is absent or is <=0, pivoting is done on all columns of A. This is appropriate if A may not be of full rank (i.e. certain columns of A are linear combinations of other columns).

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of A, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of A whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.

TOL must not be greater or equal to 1 or less or equal than 0, otherwise the numerical rank of A is determined and the calculations to determine the condition number are not performed. If TOL is absent, the numerical rank of A is determined.

MIN_NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

MIN_NORM=true, the minimun 2-norm solutions are computed.

MIN_NORM=false or if MIN_NORM is absent, solutions are computed such that if the j-th column of A is omitted from the basis, X[j,:] is set to zero.

Further Details¶

The routine first computes a QR factorization with (partial) column pivoting on option (see below):

A * P = Q * R

, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.

R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, arank, is the effective rank of A.

This leads to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.

If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

A * P = Q * T * Z

, where Z is a n-by-n orthogonal matrix and T has the form:

[ T11 T12 ]

[ T21 T22 ]

, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a arank-by-arank upper triangular matrix.

The minimum 2-norm solution is then

X = [ P * Z’ ](:,:arank) * [ inv(T11) * Q1’ * B ]

where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first arank columns of Q.

If MIN_NORM is absent or has the value false, a solution is computed as

X = P(:,:arank) * [ inv(R11) * Q1’ * B ]

where inv(R11) is the inverse of R11 and Q1 consists of the first arank columns of Q. In this case, if the j-th column of A is omitted from the basis, X[j] is set to zero.

On input, if KRANK is present and KRANK=k, the first k columns of A are to be forced into the basis. Pivoting is performed on the last n-k columns of A.

When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when A is known to be full rank.

If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of A.

TOL is an optional argument such that 0<TOL<1. If TOL is not specified, or is outside ]0,1[, the calculations to determine the condition number of A are not performed and crude tests on R(j,j) are performed to determine the numerical rank of A. If TOL is present and is in ]0,1[, the calculations to determine the condition number are performed.

If it is possible that A may not be full rank (i.e., certain columns of A are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of A so that all elements are about the same order of magnitude.

On exit, if A or B are empty, the function returns a n-by-nb matrix filled with nan() value.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`function solve_llsq ( a, b )`¶

Purpose¶

SOLVE_LLSQ computes a solution to a real linear least squares problem:

Minimize 2-norm|| B - A * X ||

A is a m-element vector, B is a m-element right hand side vector and X is a real scalar.

The function returns the solution scalar X.

A and B are not overwritten by SOLVE_LLSQ.

Arguments¶

A (INPUT) real(stnd), dimension(:)

On entry, the m-element coefficient vector A.

B (INPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

The shape of B must verify:

size( B ) = size( A ) = m .

Further Details¶

The routine first generates a real elementary reflector H of order m, such that

H * A = D , with H’ * H = I and D’ = ( d 0 )

where d is a scalar. H is represented in the form

H = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector.

The solution X is then computed as

X = [ H * B ](1) / d

On exit, if A or B are empty, the function returns a nan() value.

`function solve_llsq ( a, b )`¶

Purpose¶

SOLVE_LLSQ computes solutions to real linear least squares problems of the form:

Minimize 2-norm|| B - A * X ||

A is a m-element vector and several right hand side vectors b can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B.

The function returns the nb-element solution vector X.

A and B are not overwritten by SOLVE_LLSQ.

Arguments¶

A (INPUT) real(stnd), dimension(:)

On entry, the m-element coefficient vector A.

B (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B.

The shape of B must verify:

size( B, 1 ) = size( A ) = m .

Further Details¶

The routine first generates a real elementary reflector H of order m, such that

H * A = D , with H’ * H = I and D’ = ( d 0 )

where d is a scalar. H is represented in the form

H = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector.

The solution vector X is then computed as

X(:) = [ H * B ](1,:) / d

On exit, if A or B are empty, the function returns a nb-vector filled with nan() value.

`subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )`¶

Purpose¶

LLSQ_QR_SOLVE computes a solution to a real linear least squares problem:

Minimize 2-norm|| B - MAT * X ||

using an orthogonal factorization with columns pivoting of MAT. MAT is a m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted. Here, B is a m-element right hand side vector and X is a n-element solution vector.

MAT and B are not overwritten by LLSQ_QR_SOLVE.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-n coefficient matrix MAT.

B (INPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

The shape of B must verify:

size( B ) = size( MAT, 1 ) = m .

X (OUTPUT) real(stnd), dimension(:)

On exit, the n-element solution vector X.

The shape of X must verify:

size( X ) = size( MAT, 2 ) = n .

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution vector X.

RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the m-element residual vector for the solution vector X, RESID = B - MAT * X.

The shape of RESID must verify:

size( RESID ) = size( B ) = size( MAT, 1 ) = m .

KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)

On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT.

When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).

On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.

TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.

On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.

If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.

MIN_NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

MIN_NORM=true, the minimun 2-norm solution is computed.

MIN_NORM=false or if MIN_NORM is absent, a solution is computed such that if the j-th column of A is omitted from the basis, X[j] is set to zero.

Further Details¶

The routine first computes a QR factorization with (partial) column pivoting on option (see below):

MAT * P = Q * R

, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.

R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.

This leads to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.

If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

MAT * P = Q * T * Z

, where Z is a n-by-n orthogonal matrix and T has the form:

[ T11 T12 ]

[ T21 T22 ]

, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.

The minimum 2-norm solution is then

X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]

where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.

If MIN_NORM is absent or has the value false, a solution is computed as

X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]

where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero.

In both cases:

The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .

The m-element residual vector, B - MAT * X, can be obtained through the optional argument RESID.

On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.

When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.

On output, if KRANK is present, it contains the effective rank of MAT, i.e., the order of the submatrix R11.

TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed and its reciprocal is return in TOL.

If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )`¶

Purpose¶

LLSQ_QR_SOLVE computes solutions to real linear least squares problems of the form:

Minimize 2-norm|| B - MAT * X ||

using an orthogonal factorization with columns pivoting of MAT. MAT is a m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.

MAT and B are not overwritten by LLSQ_QR_SOLVE.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-n coefficient matrix MAT.

B (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B.

The shape of B must verify:

size( B, 1 ) = size( MAT, 1 ) = m

size( B, 2 ) = size( X, 2 ) = nb .

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the n-by-nb solution matrix X.

The shape of X must verify:

size( X, 1 ) = size( MAT, 2 ) = n

size( X, 2 ) = size( B, 2 ) = nb .

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.

The size of RNORM must verify:

size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .

RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, the residual vectors for the solutions stored columnwise in the matrix X, RESID = B - MAT * X.

The shape of RESID must verify:

size( RESID, 1 ) = size( B, 1 ) = size( MAT, 1 ) = m

size( RESID, 2 ) = size( B, 2 ) = size( X, 2 ) = nb .

KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)

On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT.

When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).

On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.

TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.

On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.

If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.

MIN_NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

MIN_NORM=true, the minimun 2-norm solutions are computed.

MIN_NORM=false or if MIN_NORM is absent, solutions are computed such that if the j-th column of A is omitted from the basis, X[j,:] is set to zero.

Further Details¶

The routine first computes a QR factorization with (partial) column pivoting on option (see below):

MAT * P = Q * R

, here P is n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.

R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.

This leads to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.

If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

MAT * P = Q * T * Z

, where Z is a n-by-n orthogonal matrix and T has the form:

[ T11 T12 ]

[ T21 T22 ]

, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.

The minimum 2-norm solution is then

X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]

where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.

If MIN_NORM is absent or has the value false, a solution is computed as

X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]

where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero.

In both cases:

The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.

The residual matrix, B - MAT * X, can be obtained through the optional argument RESID.

On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.

When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.

On output, if KRANK is present, it contains the effective rank of MAT, i.e., the order of the submatrix R11.

TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed and its reciprocal is return in TOL.

If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )`¶

Purpose¶

LLSQ_QR_SOLVE computes a solution to a real linear least squares problem:

Minimize 2-norm|| B - VEC * X ||

VEC is a m-element vector, B is a m-element right hand side vector and X is a real scalar.

VEC and B are not overwritten by LLSQ_QR_SOLVE.

Arguments¶

VEC (INPUT) real(stnd), dimension(:)

On entry, the m-element coefficient vector VEC.

B (INPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

The shape of B must verify:

size( B ) = size( VEC ) = m .

X (OUTPUT) real(stnd)

On exit, the real solution X.

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution scalar X.

RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the m-element residual vector for the solution X, RESID = B - MAT * X.

The shape of RESID must verify:

size( RESID ) = size( B ) = size( VEC ) = m .

Further Details¶

The routine first generates a real elementary reflector H of order m, such that

H * VEC = D , with H’ * H = I and D’ = ( d 0 )

where d is a scalar. H is represented in the form

H = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector.

The solution X is then computed as

X = [ H * B ](1) / d

The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM as

2-norm|| [ H * B ](2:) ||

The residual vector, B - VEC * X, can be obtained through the optional argument RESID.

`subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )`¶

Purpose¶

LLSQ_QR_SOLVE computes solutions to real linear least squares problems of the form:

Minimize 2-norm|| B - VEC * X ||

VEC is a m-element vector and several right hand side vectors b and solution scalars x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the nb-element solution vector X, respectively.

VEC and B are not overwritten by LLSQ_QR_SOLVE.

Arguments¶

VEC (INPUT) real(stnd), dimension(:)

On entry, the m-element coefficient vector VEC.

B (INPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B.

The shape of B must verify:

size( B, 1 ) = size( VEC ) = m

size( B, 2 ) = size( X ) = nb .

X (OUTPUT) real(stnd), dimension(:)

On exit, the nb-element solution vector X.

The shape of X must verify: size( X ) = size( B, 2 ) = nb .

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norms of the residual vectors for the solutions stored in the vector X.

The size of RNORM must verify:

size( RNORM ) = size( X ) = size( B, 2 ) = nb .

RESID (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, the residual vectors for the solutions stored in the vector X, RESID = B - VEC * X.

The shape of RESID must verify:

size( RESID, 1 ) = size( B, 1 ) = size( VEC ) = m

size( RESID, 2 ) = size( B, 2 ) = size( X ) = nb .

Further Details¶

The routine first generates a real elementary reflector H of order m, such that

H * A = D , with H’ * H = I and D’ = ( d 0 )

where d is a scalar. H is represented in the form

H = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector.

The solution vector X is then computed as

X(:) = [ H * B ](1,:) / d

The 2-norm of the residual vector for the solution X[j] is given in RNORM[j] if argument RNORM is present.

The residual matrix, B - VEC * X, can be obtained through the optional argument RESID.

`subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )`¶

Purpose¶

LLSQ_QR_SOLVE2 computes a solution to a real linear least squares problem:

Minimize 2-norm|| B - MAT * X ||

using a (complete) orthogonal factorization of MAT. MAT is a m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted. Here, B is a m-element right hand side vector and X is a n-element solution vector.

MAT and B are overwritten with information generated by LLSQ_QR_SOLVE2. The (complete) orthogonal factorization of MAT is saved in arguments MAT, DIAGR, BETA, IP and TAU on output.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n coefficient matrix MAT.

On exit, MAT has been overwritten by details of its (complete) orthogonal factorization. Other parts of the factorization can be obtained if the optional arguments DIAGR, BETA, IP and TAU are present.

See Further Details.

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

On exit, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B on output.

The shape of B must verify:

size( B ) = size( MAT, 1 ) = m .

X (OUTPUT) real(stnd), dimension(:)

On exit, the n-element solution vector X.

The shape of X must verify:

size( X ) = size( MAT, 2 ) = n .

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution vector X.

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B on exit.

KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)

On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT.

When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).

On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOL must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.

TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.

On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.

If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.

MIN_NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

MIN_NORM=true, a complete orthogonal factorization of MAT and the minimun 2-norm solution is computed.

MIN_NORM=false or if MIN_NORM is absent, a QR factorization with column pivoting of MAT and a solution is computed such that if the j-th column of MAT is omitted from the basis, X[j] is set to zero.

DIAGR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if:

MIN_NORM=false or is absent, the diagonal elements of the matrix R.

MIN_NORM=true, the diagonal elements of the matrix T11. The diagonal elements of T11 are stored in DIAGR(1:KRANK).

See Further Details.

The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).

BETA (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the scalar factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must be min( size(MAT,1) , size(MAT,2) ).

IP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

See Further Details.

The size of IP must be size(MAT,2) = n.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Z in the complete orthogonal factorization of MAT if MIN_NORM is present and is equal to true. Otherwise, TAU is set to 0.

See Further Details.

The size of TAU must be min( size(MAT,1) , size(MAT,2) ).

Further Details¶

The routine first computes a QR factorization with (partial) column pivoting on option (see below):

MAT * P = Q * R

, here P is a n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(k), where k = min( size(MAT,1) , size(MAT,2) )

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

The elements above the diagonal of the array MAT contain the corresponding elements of the triangular matrix R. The elements on the diagonal of R are stored in the array DIAGR.

R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.

This leads to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.

If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

MAT * P = Q * T * Z

, where Z is a n-by-n orthogonal matrix and T is a m-by-n matrix has the form:

[ T11 T12 ]

[ T21 T22 ]

, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.

The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R, is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-KRANK) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(KRANK)

On exit, the scalar tau is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,KRANK+1:n).

On exit, the elements above the diagonal of the array section MAT(1:KRANK,1:KRANK) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array section DIAGR(1:KRANK). The last part of DIAGR is set to zero. In other words, T11 overwrites R11 and Z overwrites R12 on exit.

The minimum 2-norm solution is then

X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]

where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.

If MIN_NORM is absent or has the value false, a solution is computed as

X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]

where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero and R is not destroyed in MAT.

In both cases:

The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .

If COMP_RESID=true, The m-element residual vector B - MAT * X overwrites B on exit.

On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.

When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.

On output, if KRANK is present, it contains the rank of MAT, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT.

TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed, the effective rank of MAT is determined and the reciprocal of the condition number is returned in TOL.

If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )`¶

Purpose¶

LLSQ_QR_SOLVE2 computes solutions to real linear least squares problems of the form:

Minimize 2-norm|| B - MAT * X ||

using an orthogonal factorization with columns pivoting of MAT. MAT is a m-by-n matrix which may be rank-deficient. m>=n or n>m is permitted.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.

MAT and B are overwritten with information generated by LLSQ_QR_SOLVE2. The (complete) orthogonal factorization of MAT is saved in arguments MAT, DIAGR, BETA, IP and TAU on output.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n coefficient matrix MAT.

On exit, MAT has been overwritten by details of its (complete) orthogonal factorization. Other parts of the factorization can be obtained if the optional arguments DIAGR, BETA, IP and TAU are present.

See Further Details.

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B.

On exit, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on output.

The shape of B must verify:

size( B, 1 ) = size( MAT, 1 ) = m

size( B, 2 ) = size( X, 2 ) = nb .

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the n-by-nb solution matrix X.

The shape of X must verify:

size( X, 1 ) = size( MAT, 2 ) = n

size( X, 2 ) = size( B, 2 ) = nb .

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.

The size of RNORM must verify:

size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on exit.

KRANK (INPUT/OUTPUT, OPTIONAL) integer(i4b)

On entry, KRANK=k implies that the first k columns of MAT are to be forced into the basis, pivoting is performed on the last n-k columns of MAT.

When KRANK >=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is <=0, pivoting is done on all columns of MAT. This is appropriate if MAT may not be of full rank (i.e. certain columns of MAT are linear combinations of other columns).

On exit, KRANK contains the effective rank of MAT, i.e., the number of independent columns in matrix MAT.

TOL (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of MAT whose estimated condition number, in the 1-norm, is less than 1/TOL. TOL must be set to the relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.

TOL must not be greater or equal to 1 or less than 0, otherwise the numerical rank of MAT is determined and the calculations to determine the condition number are not performed. If TOL=0, the numerical rank of MAT is determined, but the condition number is calculated.

On exit, if a condition number is calculated, its reciprocal is returned in TOL. Otherwise, TOL is not changed.

If TOL is absent, the numerical rank of MAT is used and is returned in the optional argument KRANK.

MIN_NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

MIN_NORM=true, a complete orthogonal factorization of MAT and the minimun 2-norm solutions are computed.

MIN_NORM=false or if MIN_NORM is absent, a QR factorization with column pivoting of MAT and solutions are computed such that if the j-th column of MAT is omitted from the basis, X[j,:] is set to zero.

DIAGR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the diagonal elements of the matrix R if MIN_NORM=false or is absent, or the diagonal elements of the matrix T11 if MIN_NORM is present and is equal to true. The diagonal elements of T11 are stored in DIAGR(1:KRANK).

See Further Details.

The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).

BETA (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Q.

See Further Details.

The size of BETA must be min( size(MAT,1) , size(MAT,2) ).

IP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On exit, if IP(j)=k, then the j-th column of MAT * P was the k-th column of MAT.

See Further Details.

The size of IP must be size(MAT,2) = n.

TAU (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the scalars factors of the elementary reflectors defining Z in the complete orthogonal factorization of MAT if MIN_NORM is present and is equal to true. Otherwise, TAU is set to 0.

See Further Details.

The size of TAU must be min( size(MAT,1) , size(MAT,2) ).

Further Details¶

The routine first computes a QR factorization with (partial) column pivoting on option (see below):

MAT * P = Q * R

, here P is a n-by-n permutation matrix, R is an upper triangular or trapezoidal (if n>m) matrix and Q is a m-by-m orthogonal matrix.

The matrix Q is represented as a product of elementary reflectors

Q = H(1) * H(2) * … * H(k), where k = min( size(MAT,1) , size(MAT,2) )

Each H(i) has the form

H(i) = I + beta * ( v * v’ ) ,

where beta is a real scalar and v is a real m-element vector with v(1:i-1) = 0. v(i:m) is stored on exit in MAT(i:m,i) and beta in BETA(i).

The matrix P is represented in the array IP as follows: If IP(j) = i then the jth column of P is the ith canonical unit vector.

The elements above the diagonal of the array MAT contain the corresponding elements of the triangular matrix R. The elements on the diagonal of R are stored in the array DIAGR.

R can then be partioned by defining R11 as the largest leading submatrix of R whose estimated condition number, in the 1-norm, is less than 1/TOL (or such that abs(R[j,j])>0 if TOL is absent). The order of R11, KRANK, is the effective rank of MAT.

This leads to the following partition of R:

[ R11 R12 ]

[ R21 R22 ]

where R21 is zero by construction (since R is an upper triangular or trapezoidal) and R22 is considered to be negligible.

If MIN_NORM is present and has the value true, R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

MAT * P = Q * T * Z

, where Z is a n-by-n orthogonal matrix and T has the form:

[ T11 T12 ]

[ T21 T22 ]

, here T21 (=R21), T12 and T22 (=R22) are zero and T11 is a KRANK-by-KRANK upper triangular matrix.

The factorization is obtained by Householder’s method. The kth transformation matrix, Z(k), which is used to introduce zeros into the kth row of R, is given in the form

[ I 0 ]

[ 0 T(k) ]

where

T(k) = I + tau * ( u(k) * u(k)’ ) and u(k)’ = ( 1 0 z(k) )

tau is a scalar, u(k) is a n-k+1 vector and z(k) is an (n-KRANK) element vector. tau and z(k) are chosen to annihilate the elements of the kth row of R12.

The Z n-by-n orthogonal matrix is given by

Z = Z(1) * Z(2) * … * Z(KRANK)

On exit, the scalar tau is returned in the kth element of TAU and the vector u(k) in the kth row of MAT, such that the elements of z(k) are in MAT(k,KRANK+1:n).

On exit, the elements above the diagonal of the array section MAT(1:KRANK,1:KRANK) contain the corresponding elements of the triangular matrix T11. The elements of the diagonal of T11 are stored in the array section DIAGR(1:KRANK). The last part of DIAGR is set to zero. In other words, T11 overwrites R11 and Z overwrites R12 on exit.

The minimum 2-norm solution is then

X = [ P * Z’ ](:,:KRANK) * [ inv(T11) * Q1’ * B ]

where inv(T11) is the inverse of T11, Z’ is the transpose of Z and Q1 consists of the first KRANK columns of Q.

If MIN_NORM is absent or has the value false, a solution is computed as

X = P(:,:KRANK) * [ inv(R11) * Q1’ * B ]

where inv(R11) is the inverse of R11 and Q1 consists of the first KRANK columns of Q. In this case, if the j-th column of MAT is omitted from the basis, X[j] is set to zero and R is not destroyed in MAT.

In both cases:

The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.

If COMP_RESID=true, The residual matrix B - MAT * X overwrites B on exit.

On input, if KRANK is present and KRANK=k, the first k columns of MAT are to be forced into the basis. Pivoting is performed on the last n-k columns of MAT.

When KRANK is present and KRANK>=min(m,n) is used, pivoting is not performed. This is appropriate when MAT is known to be full rank.

If KRANK is absent or is present with KRANK<=0, pivoting is done on all columns of MAT.

On output, if KRANK is present, it contains the rank of MAT, i.e., the order of the submatrix R11. This is the same as the order of the submatrix T11 in the complete orthogonal factorization of MAT.

TOL is an optional argument such that 0<=TOL<1. If TOL is not specified, or is outside [0,1[, the calculations to determine the condition number of MAT are not performed and crude tests on R(j,j) are performed to determine the numerical rank of MAT. If TOL is present and is in [0,1[, the calculations to determine the condition number are performed, the effective rank of MAT is determined and the reciprocal of the condition number is returned in TOL.

If it is possible that MAT may not be full rank (i.e., certain columns of MAT are linear combinations of other columns), then the linearly dependent columns can usually be eliminated by using KRANK=0 and TOL=relative precision of the elements in MAT and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used. Also, it may be helpful to scale the columns of MAT so that all elements are about the same order of magnitude.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )`¶

Purpose¶

LLSQ_QR_SOLVE2 computes a solution to a real linear least squares problem:

Minimize 2-norm|| B - VEC * X ||

VEC is a m-element vector, B is a m-element right hand side vector and X is a real scalar.

VEC and B are overwritten with information generated by LLSQ_QR_SOLVE2.

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element coefficient vector VEC.

On exit, VEC contains the vector v of the Householder reflector H.

See Further Details.

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

On exit, if COMP_RESID is present and is equal true, the residual vector B - VEC * X overwrites B on output.

The shape of B must verify:

size( B ) = size( VEC ) = m .

X (OUTPUT) real(stnd)

On exit, the real solution X.

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution scalar X.

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual vector B - VEC * X overwrites B on exit.

DIAGR (OUTPUT, OPTIONAL) real(stnd)

On exit, the scalar DIAGR.

See Further Details.

BETA (OUTPUT, OPTIONAL) real(stnd)

On exit, the scalar factor BETA of the elementary reflector defining H.

See Further Details.

Further Details¶

The routine first generates a real elementary reflector H of order m, such that

H * VEC = D , with H’ * H = I and D’ = ( DIAGR 0 )

where DIAGR is scalar. H is represented in the form

H = I + BETA * ( v * v’ ) ,

where BETA is a real scalar and v is a real m-element vector.

The solution X is then computed as

X = [ H * B ](1) / DIAGR

The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM as

2-norm|| [ H * B ](2:) ||

If COMP_RESID=true, The residual vector B - VEC * X overwrites B on exit.

`subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )`¶

Purpose¶

LLSQ_QR_SOLVE2 computes solutions to real linear least squares problems of the form:

Minimize 2-norm|| B - VEC * X ||

here VEC is a m-element vector and several right hand side vectors b and solution scalars x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the nb-element solution vector X, respectively.

VEC and B are overwritten with information generated by LLSQ_QR_SOLVE2.

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element coefficient vector VEC.

On exit, VEC contains the vector v of the Householder reflector H.

See Further Details.

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B.

On exit, if COMP_RESID is present and is equal true, the residual matrix B - VEC * X overwrites B on output.

The shape of B must verify:

size( B, 1 ) = size( VEC ) = m

size( B, 2 ) = size( X ) = nb .

X (OUTPUT) real(stnd), dimension(:)

On exit, the nb-element solution vector X.

The shape of X must verify:

size( X ) = size( B, 2 ) = nb .

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norms of the residual vectors for the solutions stored in the vector X.

The size of RNORM must verify:

size( RNORM ) = size( X ) = size( B, 2 ) = nb .

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual matrix B - VEC * X overwrites B on exit.

DIAGR (OUTPUT, OPTIONAL) real(stnd)

On exit, the scalar DIAGR.

See Further Details.

BETA (OUTPUT, OPTIONAL) real(stnd)

On exit, the scalar factor BETA of the elementary reflector defining H.

See Further Details.

Further Details¶

The routine first generates a real elementary reflector H of order m, such that

H * VEC = D , with H’ * H = I and D’ = ( DIAGR 0 )

where DIAGR is scalar. H is represented in the form

H = I + BETA * ( v * v’ ) ,

where BETA is a real scalar and v is a real m-element vector.

The solution vector X is then computed as

X(:) = [ H * B ](1,:) / DIAGR

The 2-norm of the residual vector for the solution X[j] is given in RNORM[j] if argument RNORM is present.

If COMP_RESID=true, The residual matrix B - VEC * X overwrites B on exit.

`subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )`¶

Purpose¶

QR_SOLVE solves overdetermined or underdetermined real linear systems

MAT * X = B

with a m-by-n matrix MAT, using a QR factorization of MAT as computed by QR_CMP. m>=n or n>m is permitted, but it is assumed that MAT has full rank. B is a m-element right hand side vector and X is a n-element solution vector.

It is assumed that QR_CMP has been used to compute the QR factorization of MAT before QR_SOLVE.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the QR factorization of the real coefficient matrix MAT as returned by QR_CMP.

The elements above the diagonal of the array contain the corresponding elements of R. The elements on and below the diagonal, with the array BETA, represent the orthogonal matrix Q in the QR decomposition of MAT, as a product of elementary reflectors, as returned by QR_CMP.

DIAGR (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the matrix R in the QR decomposition of MAT.

The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).

BETA (INPUT) real(stnd), dimension(:)

On entry, the scalars factors of the elementary reflectors defining Q, as returned by QR_CMP.

The size of BETA must be min( size(MAT,1) , size(MAT,2) ),

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

On exit, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B.

The size of B must verify:

size( B ) = size( MAT, 1 ) = m .

X (OUTPUT) real(stnd), dimension(:)

On exit, the n-element solution vector X.

The size of X must verify:

size( X ) = size( MAT, 2 ) = n .

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution vector X.

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual vector B - MAT * X overwrites B on exit.

Further Details¶

It is assumed that QR_CMP has been used to compute the QR factorization of MAT before calling QR_SOLVE.

If m>=n: the subroutine finds the least squares solution of an overdetermined system, i.e., solves the least squares problem

Minimize 2-norm|| B - MAT * X ||

If m<n: the subroutine finds a solution of an underdetermined system

MAT * X = B

The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .

If COMP_RESID=true, The m-element residual vector B - MAT * X overwrites B on exit.

MAT, DIAGR, BETA are not modified by this routine and can be left in place for successive calls with different right-hand side vectors B.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )`¶

Purpose¶

QR_SOLVE solves overdetermined or underdetermined real linear systems of the form:

MAT * X = B

with a m-by-n matrix MAT, using a QR factorization of MAT as computed by QR_CMP. m>=n or n>m is permitted, but it is assumed that MAT has full rank.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.

It is assumed that QR_CMP has been used to compute the QR factorization of MAT before QR_SOLVE.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the QR factorization of the real coefficient matrix MAT as returned by QR_CMP.

The elements above the diagonal of the array contain the corresponding elements of R. The elements on and below the diagonal, with the array BETA, represent the orthogonal matrix Q in the QR decomposition of MAT, as a product of elementary reflectors, as returned by QR_CMP.

DIAGR (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the matrix R in the QR decomposition of MAT.

The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).

BETA (INPUT) real(stnd), dimension(:)

On entry, the scalars factors of the elementary reflectors defining Q, as returned by QR_CMP.

The size of BETA must be min( size(MAT,1) , size(MAT,2) ),

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B.

On exit, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on output.

The shape of B must verify:

size( B, 1 ) = size( MAT, 1 ) = m

size( B, 2 ) = size( X, 2 ) = nb .

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the n-by-nb solution matrix X.

The shape of X must verify:

size( X, 1 ) = size( MAT, 2 ) = n

size( X, 2 ) = size( B, 2 ) = nb .

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.

The size of RNORM must verify:

size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual matrix B - MAT * X overwrites B on exit.

Further Details¶

It is assumed that QR_CMP has been used to compute the QR factorization of MAT before calling QR_SOLVE.

If m>=n: the subroutine finds the least squares solutions of overdetermined systems, i.e., solves least squares problems of the form

Minimize 2-norm|| B - MAT * X ||

If m<n: the subroutine finds solutions of underdetermined systems of the form

MAT * X = B

In both cases, several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.

The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.

If COMP_RESID=true, The residual matrix B - MAT * X overwrites B on exit.

MAT, DIAGR, BETA are not modified by this routine and can be left in place for successive calls with different right-hand side matrices B.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )`¶

Purpose¶

QR_SOLVE2 solves overdetermined or underdetermined real linear systems

MAT * X = B

with a m-by-n matrix MAT, using a QR or (complete) orthogonal factorization of MAT as computed by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines. m>=n or n>m is permitted and MAT may be rank-deficient. B is a m-element right hand side vector and X is a n-element solution vector.

It is assumed that QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 have been used to compute the (complete) orthogonal factorization of MAT before calling QR_SOLVE2.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, details of the QR or (complete) orthogonal factorization of the real coefficient matrix MAT as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

DIAGR (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the matrix R in the QR decomposition with column pivoting of MAT if TAU is absent or the diagonal elements of the matrix T11 in the complete orthogonal factorization of MAT if TAU is present, as computed by QR_CMP2. If a complete orthogonal factorization has been computed, the diagonal elements of T11 are stored in DIAGR(1:KRANK).

See description of QR_CMP2 subroutine for further details.

The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).

BETA (INPUT) real(stnd), dimension(:)

On entry, the scalars factors of the elementary reflectors defining Q in the QR or orthogonal factorization of MAT, as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

The size of BETA must be min( size(MAT,1) , size(MAT,2) ),

IP (INPUT) integer(i4b), dimension(:)

On entry, the permutation P in the QR or (complete) orthogonal factorization of MAT, as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

The size of IP must be size( MAT, 2 ) = n.

KRANK (INPUT) integer(i4b)

On entry, KRANK contains the effective rank of MAT, as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B .

On exit, if COMP_RESID is present and is equal to true, the residual vector B - MAT * X overwrites B .

The size of B must verify:

size( B ) = size( MAT, 1 ) = m .

X (OUTPUT) real(stnd), dimension(:)

On exit, the n-element solution vector X.

The size of X must verify:

size( X ) = size( MAT, 2 ) = n .

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution vector X.

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal to true, the residual vector B - MAT * X overwrites B on exit.

TAU (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a complete orthogonal factorization of MAT has been computed by QR_CMP2 and TAU contains the scalars factors of the elementary reflectors defining Z in this decomposition. Otherwise, only a QR factorization with column pivoting of MAT has been computed by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

The size of TAU must be min( size(MAT,1) , size(MAT,2) ).

Further Details¶

It is assumed that QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 has been used to compute the (complete) orthogonal factorization (if TAU is present) or the QR factorization with column pivoting (if TAU is absent) of MAT before calling QR_SOLVE2.

If m>=n: the subroutine finds the least squares solution of an overdetermined system, i.e., solves the least squares problem

Minimize 2-norm|| B - MAT * X ||

If m<n: the subroutine finds a solution of an underdetermined system

MAT * X = B

In both cases, the minimun 2-norm solution is computed if TAU is present. Otherwise, a solution is computed such that if the j-th column of MAT is omitted from the basis, X[j] is set to zero.

The 2-norm of the residual vector for the solution X can be obtained through the optional argument RNORM .

If COMP_RESID=true, The m-element residual vector B - MAT * X overwrites B on exit.

MAT, DIAGR, BETA, IP, KRANK and TAU are not modified by this routine and can be left in place for successive calls with different right-hand side vectors B.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )`¶

Purpose¶

QR_SOLVE2 solves overdetermined or underdetermined real linear systems of the form:

MAT * X = B

with a m-by-n matrix MAT, using a QR or (complete) orthogonal factorization of MAT as computed by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP and PARTIAL_RQR_CMP2 subroutines. m>=n or n>m is permitted and MAT may be rank-deficient.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.

It is assumed that QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 have been used to compute the (complete) orthogonal factorization of MAT before calling QR_SOLVE2.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, details of the QR or (complete) orthogonal factorization of the real coefficient matrix MAT as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

DIAGR (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the matrix R in the QR decomposition with column pivoting of MAT if TAU is absent or the diagonal elements of the matrix T11 in the complete orthogonal factorization of MAT if TAU is present, as computed by QR_CMP2. If a complete orthogonal factorization has been computed, the diagonal elements of T11 are stored in DIAGR(1:KRANK).

See description of QR_CMP2 subroutine for further details.

The size of DIAGR must be min( size(MAT,1) , size(MAT,2) ).

BETA (INPUT) real(stnd), dimension(:)

On entry, the scalars factors of the elementary reflectors defining Q in the QR or orthogonal factorization of MAT, as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

The size of BETA must be min( size(MAT,1) , size(MAT,2) ),

IP (INPUT) integer(i4b), dimension(:)

On entry, the permutation P in the QR or (complete) orthogonal factorization of MAT, as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

The size of IP must be size(MAT,2) = n.

KRANK (INPUT) integer(i4b)

On entry, KRANK contains the effective rank of MAT, as returned by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine for further details.

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-nb right hand side matrix B .

On exit, if COMP_RESID is present and is equal to true, the residual matrix B - MAT * X overwrites B .

The shape of B must verify:

size( B, 1 ) = size( MAT, 1 ) = m

size( B, 2 ) = size( X, 2 ) = nb .

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the n-by-nb solution matrix X.

The shape of X must verify:

size( X, 1 ) = size( MAT, 2 ) = n

size( X, 2 ) = size( B, 2 ) = nb .

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norm of the residual vectors for the solutions stored columnwise in the matrix X.

The size of RNORM must verify:

size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nb .

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal to true, the residual matrix B - MAT * X overwrites B on exit.

TAU (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a complete orthogonal factorization of MAT has been computed by QR_CMP2 and TAU contains the scalars factors of the elementary reflectors defining Z in this decomposition. Otherwise, only a QR factorization with column pivoting of MAT has been computed by QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 subroutines.

See description of QR_CMP2 subroutine further details.

The size of TAU must be min( size(MAT,1) , size(MAT,2) ).

Further Details¶

It is assumed that QR_CMP2, PARTIAL_QR_CMP, PARTIAL_RQR_CMP or PARTIAL_RQR_CMP2 has been used to compute the (complete) orthogonal factorization (if TAU is present) or the QR factorization with column pivoting (if TAU is absent) of MAT before calling QR_SOLVE2.

If m>=n: the subroutine finds the least squares solutions of overdetermined systems, i.e., solves least squares problems of the form

Minimize 2-norm|| B - MAT * X ||

If m<n: the subroutine finds solutions of underdetermined systems of the form

MAT * X = B

In both cases, several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nb right hand side matrix B and the n-by-nb solution matrix X, respectively.

In both cases, the minimun 2-norm solutions are computed if TAU is present. Otherwise, solutions are computed such that if the j-th column of MAT is omitted from the basis, X[j,:] is set to zero.

The 2-norm of the residual vector for the solution in the j-th column of X is given in RNORM[j] if argument RNORM is present.

If COMP_RESID=true, The residual matrix B - MAT * X overwrites B on exit.

MAT, DIAGR, BETA, IP, KRANK and TAU are not modified by this routine and can be left in place for successive calls with different right-hand side matrices B.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine rqb_solve ( q, b, c, x, ip, tau, comp_resid )`¶

Purpose¶

RQB_SOLVE solves overdetermined or underdetermined real linear systems

MAT * X = (Q * B) * X = C

with a m-by-n matrix MAT, using a randomized QR or complete orthogonal factorization of MAT as computed by RQB_CMP. m>=n or n>m is permitted. Q is a m-by-nqb matrix with orthonormal columns and B is a nqb-by-n upper trapezoidal matrix as computed by RQB_CMP with argument COMP_QR equals to true or with optional arguments IP and/or TAU present. C is a m-element right hand side vector and X is an approximate n-element solution vector.

It is assumed that RQB_CMP has been used to compute the randomized QR or complete orthogonal factorization of MAT before calling RQB_SOLVE.

Arguments¶

Q (INPUT) real(stnd), dimension(:,:)

On entry, the computed m-by-nqb orthonormal matrix of the randomized partial QR or complete orthogonal factorization of MAT as computed by RQB_CMP with COMP_QR equals to true or with optional arguments IP and/or TAU present.

See Further Details.

The shape of Q must verify:

size( Q, 1 ) = m ,

size( Q, 2 ) = nqb .

B (INPUT) real(stnd), dimension(:,:)

On entry, the upper trapezoidal matrix R (or T) of the randomized partial QR (or complete orthogonal) factorization of MAT as computed by RQB_CMP with COMP_QR equals to true or with optional arguments IP and/or TAU present.

See Further Details.

The shape of B must verify:

size( B, 1 ) = size( Q, 2 ) = nqb ,

size( B, 2 ) = size( X ) = n .

C (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector C.

On exit, if COMP_RESID is present and is equal true, the approximate residual vector C - MAT * X overwrites C.

The size of C must verify:

size( C ) = size( Q, 1 ) = m .

X (OUTPUT) real(stnd), dimension(:)

On exit, the approximate n-element solution vector X.

The size of X must verify:

size( X ) = size( B, 2 ) = n .

IP (INPUT, OPTIONAL) integer(i4b), dimension(:)

On entry, if IP is present a randomized partial QR or complete orthogonal factorization with column pivoting of MAT has been computed by RQB_CMP (e.g., the IP argument has also been specified in the call of RQB_CMP). If IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

If IP is present, RQB_SOLVE will reorder the elements of the solution vector X to compensate for the interchanges performed in the column pivoting phase of the QR or orthogonal factorization as computed by RQB_CMP.

See Further Details.

The size of IP must verify:

size( IP ) = size( X ) = size( B, 2 ) = n.

TAU (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a randomized partial complete orthogonal factorization of MAT has been computed by RQB_CMP (e.g., the TAU argument has also been specified in the call of RQB_CMP) and TAU stores the scalars factors of the elementary reflectors defining Z in the orthogonal factorization of MAT as computed by RQB_CMP.

If TAU is present, RQB_SOLVE will compute the approximate minimal 2-norm solution vector of the linear least-squares problem with the help of the complete orthogonal factorization of MAT as computed by RQB_CMP.

See Further Details.

The size of TAU must verify:

size( TAU ) = size( Q, 2 ) = size( B, 1 ) = nqb.

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual vector C - MAT * X overwrites C on exit.

Further Details¶

It is assumed that RQB_CMP has been used to compute the randomized partial QR or complete orthogonal factorization of MAT before calling RQB_SOLVE. RQB_CMP must be called with COMP_QR argument equals to true or with optional arguments IP and TAU eventually present. IF IP or TAU have been specified in the call of RQB_CMP, they must be also specified in the call of RQB_SOLVE.

If m>=n: the subroutine finds an approximate least squares solution of an overdetermined system, i.e., solves the least squares problem

Minimize 2-norm|| C - MAT * X ||

If m<n: the subroutine finds an approximate solution of an underdetermined system

MAT * X = C

If IP is present, RQB_SOLVE will reorder the elements of the solution vector X to compensate for the interchanges performed in the column pivoting phase of the QR or orthogonal factorization as computed by RQB_CMP.

If TAU is present, RQB_SOLVE will compute the (approximate) minimal 2-norm solution of the above least squares problems with the help of the complete orthogonal factorization of MAT as computed by RQB_CMP.

If COMP_RESID=true, The m-element residual vector C - MAT * X overwrites C on exit.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine rqb_solve ( q, b, c, x, ip, tau, comp_resid )`¶

Purpose¶

RQB_SOLVE solves overdetermined or underdetermined real linear systems

MAT * X = (Q * B) * X = C

with a m-by-n matrix MAT, using a randomized QR or complete orthogonal factorization of MAT as computed by RQB_CMP. m>=n or n>m is permitted. Q is a m-by-nqb matrix with orthonormal columns and B is a nqb-by-n upper trapezoidal matrix as computed by RQB_CMP with argument COMP_QR equals to true or with optional arguments IP and/or TAU present. C is a m-by-nc right hand side matrix and X is an approximate n-by-nc solution matrix.

It is assumed that RQB_CMP has been used to compute the randomized QR or complete orthogonal factorization of MAT before calling RQB_SOLVE.

Arguments¶

Q (INPUT) real(stnd), dimension(:,:)

On entry, the computed m-by-nqb orthonormal matrix of the randomized partial QR or complete orthogonal factorization of MAT as computed by RQB_CMP with COMP_QR equals to true or with optional arguments IP and/or TAU present.

See Further Details.

The shape of Q must verify:

size( Q, 1 ) = m ,

size( Q, 2 ) = nqb .

B (INPUT) real(stnd), dimension(:,:)

On entry, the upper trapezoidal matrix R (or T) of the randomized partial QR (or complete orthogonal) factorization of MAT as computed by RQB_CMP with COMP_QR equals to true or with optional arguments IP and/or TAU present.

See Further Details.

The shape of B must verify:

size( B, 1 ) = size( Q, 2 ) = nqb ,

size( B, 2 ) = size( X, 1 ) = n .

C (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-nc right hand side matrix C.

On exit, if COMP_RESID is present and is equal true, the approximate residual matrix C - MAT * X overwrites C.

The shape of C must verify:

size( C, 1 ) = size( Q, 1 ) = m .

size( C, 2 ) = size( X, 2 ) = nc .

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the approximate n-by-nc solution matrix X.

The shape of X must verify:

size( X, 1 ) = size( B, 2 ) = n

size( X, 2 ) = size( C, 2 ) = nc .

IP (INPUT, OPTIONAL) integer(i4b), dimension(:)

On entry, if IP is present a randomized partial QR or complete orthogonal factorization with column pivoting of MAT has been computed by RQB_CMP (e.g., the IP argument has also been specified in the call of RQB_CMP). If IP(j)=k, then the j-th column of MAT*P was the k-th column of MAT.

If IP is present, RQB_SOLVE will reorder the rows of the solution matrix X to compensate for the interchanges performed in the column pivoting phase of the QR or orthogonal factorization as computed by RQB_CMP.

See Further Details.

The size of IP must verify:

size( IP ) = size( X ) = size( B, 2 ) = n.

TAU (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if TAU is present, a randomized partial complete orthogonal factorization of MAT has been computed by RQB_CMP (e.g., the TAU argument has also been specified in the call of RQB_CMP) and TAU stores the scalars factors of the elementary reflectors defining Z in the orthogonal factorization of MAT as computed by RQB_CMP.

If TAU is present, RQB_SOLVE will compute the approximate minimal 2-norm solution matrix of the linear least-squares problem with the help of the complete orthogonal factorization of MAT as computed by RQB_CMP.

See Further Details.

The size of TAU must verify:

size( TAU ) = size( Q, 2 ) = size( B, 1 ) = nqb.

COMP_RESID (INPUT, OPTIONAL) logical(lgl)

On entry, if COMP_RESID is present and is equal true, the residual matrix C - MAT * X overwrites C on exit.

Further Details¶

It is assumed that RQB_CMP has been used to compute the randomized partial QR or complete orthogonal factorization of MAT before calling RQB_SOLVE. RQB_CMP must be called with COMP_QR argument equals to true or with optional arguments IP and TAU eventually present. IF IP or TAU have been specified in the call of RQB_CMP, they must be also specified in the call of RQB_SOLVE.

If m>=n: the subroutine finds approximate least squares solutions of overdetermined systems, i.e., solves the least squares problem

Minimize 2-norm|| C - MAT * X ||

If m<n: the subroutine finds an approximate solution of an underdetermined system

MAT * X = C

If IP is present, RQB_SOLVE will reorder the rows of the solution matrix X to compensate for the interchanges performed in the column pivoting phase of the QR or orthogonal factorization as computed by RQB_CMP.

If TAU is present, RQB_SOLVE will compute the (approximate) minimal 2-norm solutions of the above least squares problems with the help of the complete orthogonal factorization of MAT as computed by RQB_CMP.

If COMP_RESID=true, The m-by-nc residual MATRIX C - MAT * X overwrites C on exit.

For further details on linear least square problems and algorithms to solve them, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )`¶

Purpose¶

LLSQ_SVD_SOLVE computes the minimum norm solution to a real linear least squares problem:

Minimize 2-norm|| B - MAT * X ||

using the singular value decomposition (SVD) of MAT. MAT is a m-by-n matrix which may be rank-deficient.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nrhs right hand side matrix B and the n-by-nrhs solution matrix X, respectively.

The effective rank of MAT, KRANK, is determined by treating as zero those singular values which are less than TOL times the largest singular value.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

On exit, MAT is destroyed. If m>=n, MAT(:n,:n) is overwritten with the right singular vectors of MAT, stored columnwise.

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-nrhs right hand side matrix B.

On exit, B is destroyed. If m>KRANK, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements KRANK+1:m in that column.

The shape of B must verify:

size( B, 1 ) = size( MAT, 1 ) = m

size( B, 2 ) = size( X, 2 ) = nrhs .

FAILURE (OUTPUT) logical(lgl)

If:

FAILURE= false : indicates successful exit

FAILURE= true : indicates that the SVD algorithm did not converge and that full accuracy was not attained in the bidiagonal SVD of an intermediate bidiagonal form BD of MAT.

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the n-by-nrhs solution matrix X.

The shape of X must verify:

size( X, 1 ) = size( MAT, 2 ) = n

size( X, 2 ) = size( B, 2 ) = nrhs .

SINGVALUES (OUTPUT, OPTIONAL) real(stnd), dimension(:)

The singular values of MAT in decreasing order. The condition number of MAT in the 2-norm is

SINGVALUES(1)/SINGVALUES(min(m,n)).

The size of SINGVALUES must verify: size( SINGVALUES ) = min(m,n) .

KRANK (OUTPUT, OPTIONAL) integer(i4b)

On exit, the effective rank of MAT, i.e., the number of singular values which are greater than TOL * SINGVALUES(1).

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norms of the residual vectors for the solutions stored columnwise in the matrix X.

The size of RNORM must verify:

size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nrhs .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT. Singular values SINGVALUES(i) <= TOL * SINGVALUES(1) are treated as zero. If TOL is absent, machine precision is used instead.

MUL_SIZE (INPUT, OPTIONAL) integer(i4b)

Internal parameter. MUL_SIZE must verify 1 <= MUL_SIZE <= max(m,n), otherwise a default value is used. MUL_SIZE can be increased or decreased to improve the performance of the algorithm used in LLSQ_SVD_SOLVE. Maximum performance will be obtained when a real matrix of size MUL_SIZE * max(m,n) and kind stnd fits in the cache of the processors.

The default value is 32.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of MAT fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.

The default is 10.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (eg QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the singular vectors in the bidiagonal SVD algorithm.

MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of min(m,n) and the integer parameter MAX_FRANCIS_STEPS_SVD specified in the module Select_Parameters.

PERFECT_SHIFT (INPUT, OPTIONAL) logical(lgl)

PERFECT_SHIFT determines if a perfect shift strategy is used in the implicit QR algorithm in order to minimize the number of QR sweeps in the bidiagonal SVD algorithm.

The default is true.

BISECT (INPUT, OPTIONAL) logical(lgl)

BISECT determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If BISECT is set to true, singular values are computed with a more accurate bisection algorithm delivering improved accuracy in the final computed SVD decomposition and solution matrix at the expense of a slightly slower execution time.

If BISECT is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

DQDS (INPUT, OPTIONAL) logical(lgl)

DQDS determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If DQDS is set to true, singular values are computed with a more accurate dqds algorithm delivering improved accuracy in the final computed SVD decomposition at the expense of a slightly slower execution time.

If DQDS is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

Further Details¶

This subroutine is adapted from the routine DGELSS in LAPACK. If OPENMP is used, the algorithm is parallelized.

For further details on using the SVD for solving a least square problem, see the references (1), (2) or (3).

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )`¶

Purpose¶

LLSQ_SVD_SOLVE computes the minimum norm solution to a real linear least squares problem:

Minimize 2-norm|| B - MAT * X ||

using the singular value decomposition (SVD) of MAT. MAT is a m-by-n matrix which may be rank-deficient, B is a m-element right hand side vector and X is a n-element solution vector.

The effective rank of MAT, KRANK, is determined by treating as zero those singular values which are less than TOL times the largest singular value.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

On exit, MAT is destroyed. If m>=n, MAT(:n,:n) is overwritten with the right singular vectors of MAT, stored columnwise.

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

On exit, B is destroyed. If m>KRANK, the residual sum-of-squares for the solution X is given by the sum of squares of elements KRANK+1:m of B .

The size of B must verify: size( B ) = size( MAT, 1 ) = m .

FAILURE (OUTPUT) logical(lgl)

If:

FAILURE = false : indicates successful exit

FAILURE = true : indicates that the SVD algorithm did not converge and that full accuracy was not attained in the bidiagonal SVD of an intermediate bidiagonal form BD of MAT.

X (OUTPUT) real(stnd), dimension(:)

On exit, the n-element solution vector X.

The size of X must verify: size( X ) = size( MAT, 2 ) = n .

SINGVALUES (OUTPUT, OPTIONAL) real(stnd), dimension(:)

The singular values of MAT in decreasing order. The condition number of MAT in the 2-norm is

SINGVALUES(1)/SINGVALUES(min(m,n)).

The size of SINGVALUES must verify: size( SINGVALUES ) = min(m,n) .

KRANK (OUTPUT, OPTIONAL) integer(i4b)

On exit, the effective rank of MAT, i.e., the number of singular values which are greater than TOL * SINGVALUES(1).

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution vector X.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT. Singular values SINGVALUES(i) <= TOL * SINGVALUES(1) are treated as zero. If TOL is absent, machine precision is used instead.

MUL_SIZE (INPUT, OPTIONAL) integer(i4b)

Internal parameter. MUL_SIZE must verify 1 <= MUL_SIZE <= max(m,n), otherwise a default value is used. MUL_SIZE can be increased or decreased to improve the performance of the algorithm used in LLSQ_SVD_SOLVE. Maximum performance will be obtained when a real matrix of size MUL_SIZE * max(m,n) and kind stnd fits in the cache of the processors.

The default value is 32.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of MAT fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.

The default value is 10.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (eg QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the singular vectors in the bidiagonal SVD algorithm.

MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of min(m,n) and the integer parameter MAX_FRANCIS_STEPS_SVD specified in the module Select_Parameters.

PERFECT_SHIFT (INPUT, OPTIONAL) logical(lgl)

PERFECT_SHIFT determines if a perfect shift strategy is used in the implicit QR algorithm in order to minimize the number of QR sweeps in the bidiagonal SVD algorithm.

The default is true.

BISECT (INPUT, OPTIONAL) logical(lgl)

BISECT determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If BISECT is set to true, singular values are computed with a more accurate bisection algorithm delivering improved accuracy in the final computed SVD decomposition and solution vector at the expense of a slightly slower execution time.

If BISECT is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

DQDS (INPUT, OPTIONAL) logical(lgl)

DQDS determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If DQDS is set to true, singular values are computed with a more accurate dqds algorithm delivering improved accuracy in the final computed SVD decomposition at the expense of a slightly slower execution time.

If DQDS is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

Further Details¶

This subroutine is adapted from the routine DGELSS in LAPACK software. If OPENMP is used, the algorithm is parallelized.

For further details on using the SVD for solving a least square problem, see the references (1), (2) or (3).

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

`subroutine llsq_svd_solve2 ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )`¶

Purpose¶

LLSQ_SVD_SOLVE2 computes the minimum norm solution to a real linear least squares problem:

Minimize 2-norm|| B - MAT * X ||

using the one-sided Ralha-Barlow bidiagonal reduction algorithm and the singular value decomposition (SVD) of MAT. MAT is a m-by-n matrix with m >= n, which may be rank-deficient.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the m-by-nrhs right hand side matrix B and the n-by-nrhs solution matrix X, respectively.

The effective rank of MAT, KRANK, is determined by treating as zero those singular values which are less than TOL times the largest singular value.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

On exit, MAT is destroyed.

The shape of MAT must verify:

size( MAT, 1 ) >= size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-nrhs right hand side matrix B.

On exit, B is destroyed. If m>KRANK, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements KRANK+1:m in that column.

The shape of B must verify:

size( B, 1 ) = size( MAT, 1 ) = m

size( B, 2 ) = size( X, 2 ) = nrhs .

FAILURE (OUTPUT) logical(lgl)

If:

FAILURE= false : indicates successful exit

FAILURE= true : indicates that the SVD algorithm did not converge and that full accuracy was not attained in the bidiagonal SVD of an intermediate bidiagonal form BD of MAT.

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the n-by-nrhs solution matrix X.

The shape of X must verify:

size( X, 1 ) = size( MAT, 2 ) = n

size( X, 2 ) = size( B, 2 ) = nrhs .

SINGVALUES (OUTPUT, OPTIONAL) real(stnd), dimension(:)

The singular values of MAT in decreasing order. The condition number of MAT in the 2-norm is

SINGVALUES(1)/SINGVALUES(min(m,n)).

The size of SINGVALUES must verify: size( SINGVALUES ) = min(m,n) .

KRANK (OUTPUT, OPTIONAL) integer(i4b)

On exit, the effective rank of MAT, i.e., the number of singular values which are greater than TOL * SINGVALUES(1).

RNORM (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the 2-norms of the residual vectors for the solutions stored columnwise in the matrix X.

The size of RNORM must verify:

size( RNORM ) = size( X, 2 ) = size( B, 2 ) = nrhs .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT. Singular values SINGVALUES(i) <= TOL * SINGVALUES(1) are treated as zero. If TOL is absent, machine precision is used instead.

MUL_SIZE (INPUT, OPTIONAL) integer(i4b)

Internal parameter. MUL_SIZE must verify 1 <= MUL_SIZE <= max(m,n), otherwise a default value is used. MUL_SIZE can be increased or decreased to improve the performance of the algorithm used in LLSQ_SVD_SOLVE. Maximum performance will be obtained when a real matrix of size MUL_SIZE * max(m,n) and kind stnd fits in the cache of the processors.

The default value is 32.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of MAT fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.

The default is 10.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (eg QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the singular vectors in the bidiagonal SVD algorithm.

MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of min(m,n) and the integer parameter MAX_FRANCIS_STEPS_SVD specified in the module Select_Parameters.

PERFECT_SHIFT (INPUT, OPTIONAL) logical(lgl)

PERFECT_SHIFT determines if a perfect shift strategy is used in the implicit QR algorithm in order to minimize the number of QR sweeps in the bidiagonal SVD algorithm.

The default is true.

BISECT (INPUT, OPTIONAL) logical(lgl)

BISECT determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If BISECT is set to true, singular values are computed with a more accurate bisection algorithm delivering improved accuracy in the final computed SVD decomposition and solution matrix at the expense of a slightly slower execution time.

If BISECT is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

DQDS (INPUT, OPTIONAL) logical(lgl)

DQDS determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If DQDS is set to true, singular values are computed with a more accurate dqds algorithm delivering improved accuracy in the final computed SVD decomposition at the expense of a slightly slower execution time.

If DQDS is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

Further Details¶

This subroutine is adapted from the routine DGELSS in LAPACK. If OPENMP is used, the algorithm is parallelized.

This subroutine is faster than the LLSQ_SVD_SOLVE subroutine for solving linear least square problems with a huge coefficient matrix thanks to the use of the Ralha-Barlow one-sided bidiagonal reduction algorithm, but is restricted to the case where the coefficient matrix MAT has more rows than columns.

For further details on the Ralha-Barlow one-sided bidiagonal reduction algorithm and using the SVD for solving a least square problem, see the references (1), (2), (3), (4) and (5).

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

Barlow, J.L., Bosner, N., and Drmac, Z., 2005:

A new stable bidiagonal reduction algorithm. Linear Algebra Appl., No 397, pp. 35-84.

Bosner, N., and Barlow, J.L., 2007:

Block and Parallel versions of one-sided bidiagonalization. SIAM J. Matrix Anal. Appl., Volume 29, No 3, pp. 927-953.

`subroutine llsq_svd_solve2 ( mat, b, failure, x, &`¶

mul_size, maxiter, max_francis_steps, & perfect_shift, bisect, dqds )

Purpose¶

LLSQ_SVD_SOLVE2 computes the minimum norm solution to a real linear least squares problem:

Minimize 2-norm|| B - MAT * X ||

using the one-sided Ralha-Barlow bidiagonal reduction algorithm and the singular value decomposition (SVD) of MAT. MAT is a m-by-n matrix with m >= n, which may be rank-deficient. B is a m-element right hand side vector and X is a n-element solution vector.

The effective rank of MAT, KRANK, is determined by treating as zero those singular values which are less than TOL times the largest singular value.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the m-by-n matrix MAT.

On exit, MAT is destroyed.

The shape of MAT must verify:

size( MAT, 1 ) >= size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the m-element right hand side vector B.

On exit, B is destroyed. If m>KRANK, the residual sum-of-squares for the solution X is given by the sum of squares of elements KRANK+1:m of B .

The size of B must verify: size( B ) = size( MAT, 1 ) = m .

FAILURE (OUTPUT) logical(lgl)

If:

FAILURE = false : indicates successful exit

FAILURE = true : indicates that the SVD algorithm did not converge and that full accuracy was not attained in the bidiagonal SVD of an intermediate bidiagonal form BD of MAT.

X (OUTPUT) real(stnd), dimension(:)

On exit, the n-element solution vector X.

The size of X must verify: size( X ) = size( MAT, 2 ) = n .

SINGVALUES (OUTPUT, OPTIONAL) real(stnd), dimension(:)

The singular values of MAT in decreasing order. The condition number of MAT in the 2-norm is

SINGVALUES(1)/SINGVALUES(min(m,n)).

The size of SINGVALUES must verify: size( SINGVALUES ) = min(m,n) .

KRANK (OUTPUT, OPTIONAL) integer(i4b)

On exit, the effective rank of MAT, i.e., the number of singular values which are greater than TOL * SINGVALUES(1).

RNORM (OUTPUT, OPTIONAL) real(stnd)

On exit, the 2-norm of the residual vector for the solution vector X.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of MAT. Singular values SINGVALUES(i) <= TOL * SINGVALUES(1) are treated as zero. If TOL is absent, machine precision is used instead.

MUL_SIZE (INPUT, OPTIONAL) integer(i4b)

Internal parameter. MUL_SIZE must verify 1 <= MUL_SIZE <= max(m,n), otherwise a default value is used. MUL_SIZE can be increased or decreased to improve the performance of the algorithm used in LLSQ_SVD_SOLVE. Maximum performance will be obtained when a real matrix of size MUL_SIZE * max(m,n) and kind stnd fits in the cache of the processors.

The default value is 32.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of MAT fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.

The default value is 10.

MAX_FRANCIS_STEPS (INPUT, OPTIONAL) integer(i4b)

MAX_FRANCIS_STEPS controls the maximum number of Francis sets (eg QR sweeps) of Givens rotations which must be saved before applying them with a wavefront algorithm to accumulate the singular vectors in the bidiagonal SVD algorithm.

MAX_FRANCIS_STEPS is a strictly positive integer, otherwise the default value is used.

The default is the minimum of min(m,n) and the integer parameter MAX_FRANCIS_STEPS_SVD specified in the module Select_Parameters.

PERFECT_SHIFT (INPUT, OPTIONAL) logical(lgl)

PERFECT_SHIFT determines if a perfect shift strategy is used in the implicit QR algorithm in order to minimize the number of QR sweeps in the bidiagonal SVD algorithm.

The default is true.

BISECT (INPUT, OPTIONAL) logical(lgl)

BISECT determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If BISECT is set to true, singular values are computed with a more accurate bisection algorithm delivering improved accuracy in the final computed SVD decomposition and solution vector at the expense of a slightly slower execution time.

If BISECT is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

DQDS (INPUT, OPTIONAL) logical(lgl)

DQDS determines how the singular values are computed if a perfect shift strategy is used in the bidiagonal SVD algorithm (e.g., if PERFECT_SHIFT is equal to TRUE). This argument has no effect if PERFECT_SHIFT is equal to false.

If DQDS is set to true, singular values are computed with a more accurate dqds algorithm delivering improved accuracy in the final computed SVD decomposition at the expense of a slightly slower execution time.

If DQDS is set to false, singular values are computed with the fast Pal-Walker-Kahan algorithm.

If both optional arguments BISECT and DQDS are specified with the value true, the bisection algorithm is used.

The default is false.

Further Details¶

This subroutine is adapted from the routine DGELSS in LAPACK software. If OPENMP is used, the algorithm is parallelized.

This subroutine is faster than the LLSQ_SVD_SOLVE subroutine for solving linear least square problems with a huge coefficient matrix thanks to the use of the Ralha-Barlow one-sided bidiagonal reduction algorithm, but is restricted to the case where the coefficient matrix MAT has more rows than columns.

For further details on the Ralha-Barlow one-sided bidiagonal reduction algorithm and using the SVD for solving a least square problem, see the references (1), (2), (3), (4) and (5).

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Lawson, C.L., and Hanson, R.J., 1974:

Solving least square problems. Prentice-Hall.

Hansen, P.C., Pereyra, V., and Scherer, G., 2012:

Least Squares Data Fitting with Applications. Johns Hopkins University Press, 328 pp.

Barlow, J.L., Bosner, N., and Drmac, Z., 2005:

A new stable bidiagonal reduction algorithm. Linear Algebra Appl., No 397, pp. 35-84.

Bosner, N., and Barlow, J.L., 2007:

Block and Parallel versions of one-sided bidiagonalization. SIAM J. Matrix Anal. Appl., Volume 29, No 3, pp. 927-953.

Module_LLSQ_Procedures¶

function solve_llsq ( a, b, krank, tol, min_norm )¶

Purpose¶

Arguments¶

Further Details¶

function solve_llsq ( a, b, krank, tol, min_norm )¶

Purpose¶

Arguments¶

Further Details¶

function solve_llsq ( a, b )¶

Purpose¶

Arguments¶

Further Details¶

function solve_llsq ( a, b )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )¶

Purpose¶

Arguments¶

Further Details¶

subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )¶

Purpose¶

Arguments¶

Further Details¶

subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )¶

Purpose¶

Arguments¶

Further Details¶

subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )¶

Purpose¶

Arguments¶

Further Details¶

subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )¶

Purpose¶

Arguments¶

Further Details¶

subroutine rqb_solve ( q, b, c, x, ip, tau, comp_resid )¶

Purpose¶

Arguments¶

Further Details¶

subroutine rqb_solve ( q, b, c, x, ip, tau, comp_resid )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )¶

Purpose¶

Arguments¶

Further Details¶

subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )¶

Purpose¶

Arguments¶

`function solve_llsq ( a, b, krank, tol, min_norm )`¶

`function solve_llsq ( a, b, krank, tol, min_norm )`¶

`function solve_llsq ( a, b )`¶

`function solve_llsq ( a, b )`¶

`subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )`¶

`subroutine llsq_qr_solve ( mat, b, x, rnorm, resid, krank, tol, min_norm )`¶

`subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )`¶

`subroutine llsq_qr_solve ( vec, b, x, rnorm, resid )`¶

`subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )`¶

`subroutine llsq_qr_solve2 ( mat, b, x, rnorm, comp_resid, krank, tol, min_norm, diagr, beta, ip, tau )`¶

`subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )`¶

`subroutine llsq_qr_solve2 ( vec, b, x, rnorm, comp_resid, diagr, beta )`¶

`subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )`¶

`subroutine qr_solve ( mat, diagr, beta, b, x, rnorm, comp_resid )`¶

`subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )`¶

`subroutine qr_solve2 ( mat, diagr, beta, ip, krank, b, x, rnorm, comp_resid, tau )`¶

`subroutine rqb_solve ( q, b, c, x, ip, tau, comp_resid )`¶

`subroutine rqb_solve ( q, b, c, x, ip, tau, comp_resid )`¶

`subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )`¶

`subroutine llsq_svd_solve ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )`¶

`subroutine llsq_svd_solve2 ( mat, b, failure, x, singvalues, krank, rnorm, tol, mul_size, maxiter, max_francis_steps, perfect_shift, bisect, dqds )`¶

`subroutine llsq_svd_solve2 ( mat, b, failure, x, &`¶