f04qaf
f04qaf
© Numerical Algorithms Group, 2002.
Purpose
F04QAF Sparse linear least-squares problem, m real equations in n unknowns
Synopsis
[x,se,itn,anorm,acond,rnorm,arnorm,xnorm,inform,ifail] = f04qaf(n,b,aprod,...
damp<,atol,btol,conlim,itnlim,msglvl,ifail>)
Description
F04QAF can be used to solve a system of linear equations
Ax=b (3.1)
where A is an n by n sparse unsymmetric matrix, or can be used to
solve linear least-squares problems, so that F04QAF minimizes the
value (rho) given by
(rho)=||r||, r=b-Ax (3.2)
where A is an m by n sparse matrix and ||r|| denotes the
2 T
Euclidean length of r so that ||r|| =r r. A damping parameter,
(lambda), may be included in the least-squares problem in which
case F04QAF minimizes the value (rho) given by
2 2 2 2
(rho) =||r|| +(lambda) ||x|| (3.3)
(lambda) is supplied as the parameter DAMP and should of course
be zero if the solution to problems (3.1) or (3.2) is required.
Minimizing (rho) in (3.3) is often called ridge regression.
F04QAF solves the problems by an algorithm based upon the Lanczos
process. The routine does not require A explicitly, but A is
specified via a user-supplied routine APROD which must perform
T
the operations (y+Ax) and (x+A y) for a given n element vector x
and m element vector y. A parameter to APROD specifies which of
the two operations is required on a given entry.
The routine also returns estimates of the standard errors of the
sample regression coefficients (x , for i=1,2,...,n) given by the
i
diagonal elements of the estimated variance-covariance matrix V.
When problem (3.2) is being solved and A is of full rank, then V
is given by
2 T -1 2 2
V=s (A A) , s =(rho) /(m-n), m>n
and when problem (3.3) is being solved then V is given by
2 T 2 -1 2 2
V=s (A A+(lambda) I) , s =(rho) /m, (lambda)/=0.
_
Let A denote the matrix
_ _ (A )
A=A, (lambda)=0 ; A=((lambda)I), (lambda)/=0, (3.4)
_
let r denote the residual vector
_ _ (b) _
r=r, (lambda)=0 ; r=(0)-Ax, (lambda)/=0 (3.5)
_
corresponding to an iterate x, so that (rho)=||r|| is the
function being minimized, and let ||A|| denote the Frobenius
(Euclidean) norm of A. Then the routine accepts x as a solution
if it is estimated that one of the following two conditions is
satisfied:
_
(rho)<=tol ||A||.||x||+tol ||b|| (3.6)
1 2
_T_ _
||A r||<=tol ||A||(rho) (3.7)
1
where tol and tol are user-supplied tolerances which estimate
1 2
the relative errors in A and b respectively. Condition (3.6) is
appropriate for compatible problems where, in theory, we expect
the residual to be zero and will be satisfied by an acceptable
solution x to a compatible problem. Condition (3.7) is
appropriate for incompatible systems where we do not expect the
residual to be zero and is based upon the observation that, in
theory,
_T_
A r=0
when x is a solution to the least-squares problem, and so (3.7)
will be satisfied by an acceptable solution x to a linear least-
squares problem.
The routine also includes a test to prevent convergence to
solutions, x, with unacceptably large elements. This can happen
if A is nearly singular or is nearly rank deficient. If we let
_
the singular values of A be
(sigma) >=(sigma) >=...>=(sigma) >=0
1 2 n
_
then the condition number of A is defined as
_
cond(A)=(sigma) /(sigma)
1 k
_
where (sigma) is the smallest non-zero singular value of A and
k
_ _
hence k is the rank of A. When k<n, then A is rank deficient, the
least-squares solution is not unique and F04QAF will normally
_
converge to the minimal length solution. In practice A will not
have exactly zero singular values, but may instead have small
singular values that we wish to regard as zero.
The routine provides for this possibility by terminating if
_
cond(A)>=c (3.8)
lim
_
where c is a user-supplied limit on the condition number of A.
lim
For problem (3.1) termination with this condition indicates that
A is nearly singular and for problem (3.2) indicates that A is
nearly rank deficient and so has near linear dependencies in its
T
columns. In this case inspection of ||r||, ||A r|| and ||x||,
which are all returned by the routine, will indicate whether or
not an acceptable solution has been found. Condition (3.8),
perhaps in conjunction with (lambda)/=0, can be used to try and '
regularise' least-squares solutions.
Introduction of a non-zero damping parameter (lambda) tends to
reduce the size of the computed solution and to make its
components less sensitive to changes in the data, and F04QAF is
applicable when a value of (lambda) is known a priori. To have an
_________
effect, (lambda) should normally be at least \/(epsilon)||A||
where (epsilon) is the machine precision.
Whenever possible the matrix A should be scaled so that the
relative errors in the elements of A are all of comparable size.
Such a scaling helps to prevent the least-squares problem from
being unnecessarily sensitive to data errors and will normally
reduce the number of iterations required. At the very least, in
the absence of better information, the columns of A should be
scaled to have roughly equal column length.
Parameters
f04qaf
Required Input Arguments:
n integer
b (:) real
aprod function (User-Supplied)
damp real
Optional Input Arguments: <Default>
atol real eps
btol real eps
conlim real 1/atol
itnlim integer 2*n
msglvl integer 0
ifail integer -1
Output Arguments:
x (n) real
se (n) real
itn integer
anorm real
acond real
rnorm real
arnorm real
xnorm real
inform integer
ifail integer