Methods for Large Scale Total Least Squares Problems.
The solution of the total least squares (TLS) problems, $\min_{E,f}\|(E,f)\|_F$ subject to (A+E)x=b+f, can in the generic case be obtained from the right singular vector corresponding to the smallest singular value $\sigma_{n+1}$ of (A, b). When A is large and sparse (or structured) a method based on Rayleigh quotient iteration (RQI) has been suggested by Bjrck. In this method the problem is reduced to the solution of a sequence of symmetric, positive definite linear systems of the form $(A^TA-\bar\sigma^2I)z=g$, where $\bar\sigma$ is an approximation to $\sigma_{n+1}$. These linear systems are then solved by a {\em preconditioned} conjugate gradient method (PCGTLS). For TLS problems where A is large and sparse a (possibly incomplete) Cholesky factor of ATA can usually be computed, and this provides a very efficient preconditioner. The resulting method can be used to solve a much wider range of problems than it is possible to solve by using Lanczos-type algorithms directly for the singular value problem. In this paper the RQI-PCGTLS method is further developed, and the choice of initial approximation and termination criteria are discussed. Numerical results confirm that the given algorithm achieves rapid convergence and good accuracy.}
Introduction
.
The estimation of parameters in linear models is a fundamental problem in
many scientific and engineering applications. A statistical model that is often
realistic is to assume that the parameters x to be determined satisfy a linear
relation
where A 2 R m\Thetan , and b 2 R m , are known and (E; f) is an error matrix with
rows which are independently and identically distributed with zero mean and
the same variance. (To satisfy this assumption the data (A; b) may need to be
premultiplied by appropriate scaling matrices, see Golub and Van Loan [10].)
In statistics this model is known as the "errors-in-variables model".
The estimate of the true but unknown parameter vector x in the model (1.1)
is obtained from the solution of the total least squares (TLS) problem
min
subject to
Department of Mathematics, University of Link-oping, S-581 83 Link-oping, Sweden. e-mail:
akbjo@math.liu.se, pontus.matstoms@vti.se. The work of these authors was supported by the
Swedish Research Council for Engineering Sciences, TFR.
y Department of Informatics, University of Bergen, N-5020 Bergen, Norway, email:
pinar@ii.uib.no
denotes the Frobenius matrix norm. If a minimizing pair (E; f)
has been found for the problem (1.2) then any x satisfying
said to solve the TLS problem.
Due to recent advances in data collection techniques LS or TLS problems
where A is large and sparse (or structured) frequently arise, e.g., in signal and
image processing applications. For the solution of the LS problem both direct
methods based on sparse matrix factorizations and iterative methods are well
developed, see [2].
An excellent treatment of theoretical and computational aspects of the TLS
problem is given in Van Huffel and Vandewalle [25]. Solving the TLS problem
requires the computation of the smallest singular value and the corresponding
right singular vector of (A; b). When A is large and sparse this is a much more
difficult problem than that of computing the LS solution. For example, it is
usually not feasible to compute the SVD or any other two-sided orthogonal
factorization of A since the factors typically are not sparse.
Iterative algorithms for computing the singular subspace of a matrix associated
with its smallest singular values, with applications to TLS problems with slowly
varying data, have previously been studied by Van Huffel [24]. In [27, 3] a new
class of methods based on a Rayleigh quotient iteration was developed for the
efficient solution of large scale TLS problems. Related methods for Toeplitz
systems were studied by Kamm and Nagy [14]. In this paper the methods in [3]
are further developed and numerical results given. Similar algorithms for solving
large scale multidimensional TLS problems will be considered in a forthcoming
paper [4].
In Section 2 we recall how the solution to the TLS problem can be expressed
in terms of the smallest singular value and corresponding right singular vector
of the compound matrix (A; b). We discuss the conditioning of the LS and TLS
problems and illustrate how the TLS problem can rapidly become intractable.
Section 3 first reviews a Newton iteration for solving a secular equation. For
this method to converge to the TLS solution strict conditions on the initial approximation
have to be satisfied. We then derive the Rayleigh quotient method,
which ultimately achieves cubic convergence. The choice of initial estimates and
termination criteria are discussed. A preconditioned conjugate gradient method
is developed in Section 4 for the efficient solution of the resulting sequence of
sparse symmetric linear systems. Finally, in Section 5, numerical results are
given which confirm the rapid convergence and numerical stability of this class
of methods.
Preliminaries.
2.1 The TLS problem.
The TLS problem (1.2) is equivalent to finding a perturbation matrix (E; f)
having minimal Frobenius norm, which lowers the rank of the matrix (A; b).
Hence it can be analyzed in terms of the singular value decomposition
are the singular values of (A; b). Note that
by the minmax characterization of singular values it follows that the singular
values oe 0
i of A interlace those of (A; b), i.e.,
We assume in the following that A has full rank, that is, oe that
. Then the minimum is attained for the rank one perturbation
for which k(E; f)k solution is then obtained from the right
singular vector
z
x TLS
provided that i 6= 0. If the TLS problem is called nongeneric, and there
is no solution. This case cannot occur if oe and in the following we
always assume that this condition holds.
From the characterization (2.2) it follows that
n+1 and
the system of nonlinear equations
A T A A T b
x
x
Putting
n+1 the first block row of this system of equations can be written
which can be viewed as "the normal equations" for the TLS problem. Note that
from our assumption that oe 0 n ? oe n+1 it follows that A T A \Gamma oe 2
n+1 I is positive
definite.
2.2 Conditioning of the TLS problem.
For the evaluation of accuracy and stability of the algorithms to be presented
we need to know the sensitivity of the TLS problem to perturbations in data.
We first recall that if x LS 6= 0 the condition number for the LS problem is (see
[2, Sec. 1.4])
n . Note that the condition number depends on both A and
b, and that for large residual problems the second term may dominate.
Condition number for problem LS and TLS
beta
kappa
kappa LS
kappa TLS

Figure

2.1: Condition numbers -LS and - TLS as function of
Equation (2.4) shows that the TLS problem is always worse conditioned than
the LS problem. From (2.3), multiplying from the left with
This inequality is wek, but shows that kx TLS k 2 will be large when kr LS k 2 AE oe 0
n .
Golub and Van Loan [10] showed that an approximate condition number for
the TLS problem is
the TLS condition number can be much greater than
-(A). The relation between the two condition numbers (2.5) and (2.7) depend
on the relation between the kr LS k 2 and oe n+1 , which is quite intricate. (For a
study of this relation in another context see Paige and Strako-s [17].)
As an illustration we consider the following small overdetermined system@
Trivially, the LS solution is
If we take in (2.8) oe independent
of fi, and hence does not reflect the illconditioning of A. The TLS solution is of
similar size as the LS solution as long as jfij - oe 0 2 . However, when jfij AE oe 0 2
then from (2.6) it follows that kx TLS k 2 is large.
In Fig. 2.1 the two condition numbers are plotted as a function of jfij. We note
that - LS increases proportionally to jfij because of the second term in (2.5). For
the condition number - TLS grows proportionally to jfij 2 . It can be
verified that kx TLS k 2 also grows proportionally to jfij 2 .
3 Newton and Rayleigh Quotient methods.
3.1 A Newton method.
Equation (2.3) constitutes a system of (n equations in x and
-. One way to proceed (see [14]) is to eliminate x to obtain the rational secular
equation for
method applied to (3.1) leads to the
iteration
This iteration will converge monotonically at a rate that is asymptotically quad-
ratic. The convergence of this method can be improved by using a rational
interpolation similar to that in [6] to solve the secular equation. However, in
any case, - will converge to oe 2
n+1 and x (k) to the TLS solution only if the initial
approximation satisfies
In general it is hard to verify this assumption. For the special case of a Toeplitz
TLS problem Kamm and Nagy [14] use a bisection algorithm based on a fast
algorithm for factorizing Toeplitz matrices to find an initial starting value satisfying
(3.4).
3.2 The Rayleigh quotient method.
The main drawback of the Newton method above is that unless (3.4) is satisfied
it will converge to the wrong singular value. A different Newton method is
obtained by applying Newton's method to the full system
As remarked in [20] this is closely related to inverse itera-
tion, which is one of the most widely used methods for refining eigenvalues and
eigenvectors. Rayleigh quotient iteration (RQI) is inverse iteration with a shift
equal to the Rayleigh quotient. RQI has cubic convergence for the symmetric
eigenvalue problem, see [18, Sec.4-7], and is superior to the standard Newton
method applied to (3.5).
For the eigenvalue problem (2.3) the Rayleigh quotient equals
Let x (k) be the current approximation and ae k the corresponding Rayleigh
quotient. Then the next approximation x (k+1) in RQI and the scaling factor fi k
are obtained from the symmetric linear system
where
If J (k) is positive definite the solution can be obtained by block Gaussian elimi-
nation, '
\Gamma(z
where
It follows that x
In [2] a reformulation was made to express the solution in terms of the residual
vectors of (3.5) '
where r This uses the following formulas to compute
The RQI iteration is defined by equations (3.10)-(3.13).
3.3 Initial estimate and global convergence.
Parlett and Kahan [19] have shown that for almost all initial vectors the
Rayleigh quotient iteration converges to some singular value and vector pair.
However, in general we cannot say to which singular vector RQI will converge.
If the LS solution is known, a suitable starting approximation for - may be
Conditions to ensure that RQI will converge to the TLS solution from the
starting approximation (ae(x LS ); x LS ) are in general difficult to verify and often
not satisfied in practice. However, in contrast to the simple Newton iteration
in Section 3.1, the method may converge to the TLS solution even when
The Rayleigh quotient ae(x LS ) will be a large overestimate of oe 2
n+1 when the
residual norm kr LS k 2 is large and kx LS k 2 does not reflect the illconditioning of
A. Note that it is typical for illconditioned least squares problems that the right-hand
side is such that kx LS k 2 is not large! For example, least squares problems
arising from ill-posed problems usually satisfy a so called Picard condition, which
guarantees that the right-hand side has this property, see [11, Sec. 1.2.3].
Szyld [23] suggested that one or more steps of inverse iteration could be applied
initially before switching to RQI, in order to ensure convergence to the smallest
eigenvalue. Inverse iteration for oe 2
n+1 corresponds to taking oe in the RQI
algorithm. Starting from x = x LS the first step of inverse iteration simplifies as
follows. Using (3.9) and (3.10) with ae
z
and the new approximation becomes
Several steps of inverse iteration may be needed to ensure convergence of RQI to
the smallest singular value. However, since inverse iteration only converges lin-
early, taking more than one step will usually just hold up the rapid convergence
of RQI. We therefore recommend in general steps as the default value.
To illustrate the situation consider again the small 3 \Theta 2 system (2.8) with
. This has the LS solution x
does not reflect the illconditioning of A the initial
Rayleigh quotient approximation equals
By the interlacing property we have that oe 3 - oe 0 2 . Since jfij AE oe 0 2 it is clear that
the Rayleigh quotient fails to approximate oe 2
3 . This is illustrated in Figure 3.1,
where ae(x LS ) 1=2 and oe 3 are plotted as function of jfij. It is easily verified,
however, that after one step of inverse iteration ae(x INV ) will be close to oe 0 2
-5
beta
sqrt(r LS

Figure

3.1: Rayleigh quotient approximation and oe 3 for
3.4 Termination criteria for RQI.
The RQI algorithm for the TLS problem is defined by (3.10)-(3.13). When
should the RQI iteration be terminated? We suggest two different criteria.
The first is based on the key fact in the proof of global convergence that the
normalized residual norm
always decreases, fl k+1 - fl k , for all k. Thus, if an increase in the norm occurs
this must be caused by roundoff, and then it makes no sense to continue the
iterations. This suggests that we terminate the iterations with x k+1 when
A second criterion is based on the observation that since the condition number
for computing oe n+1 equals 1, we can expect to obtain oe n+1 to full machine
precision. Since convergence of RQI is cubic a criterion could be to stop when
the change in the approximation to oe n+1 is of the order of oe 1 u 1=p , where
similar criterion with used by Kamm and Nagy [14] for terminating
the Newton iteration.) However, as will be evident from the numerical results
in Section 5, full accuracy in x TLS in general requires one more iteration after
oe n+1 has converged. Therefore we recommend to stop when either (3.16) or
is satisfies, where u is the machine unit and C a suitable constant.
We summarize below the RQI algorithm with one step of inverse iteration (cf.
Algorithm 3.1. Rayleigh Quotient Iteration.
solve A T
solve
solve
3.5 Rounding errors and stability.
If the RQI iteration converges then f (k) , g (k) , and fi k will tend to zero. Consider
the rounding errors which occur in the evaluation of the residuals (3.11). Let
~
where u is the unit roundoff; see [13, Chap. 3]. Then the computed
residual vector satisfies -
Obviously convergence will cease when the residuals (3.11) are dominated by
roundoff. Assume that we perform one iteration from the exact solution, x TLS ,
r TLS , and
n+1 . Then the first correction to the current approximation is
obtained by solving the linear system in (3.13), which now becomes
For the correction this gives the estimate
This estimate is consistent with the condition estimate for the TLS problem.
We note that the equations (3.18) are of similar form to those that appear in
the corrected semi-normal equations for the LS problem; see [1], [2, Sec. 6.6.5].
A detailed roundoff error analysis similar to that done for the LS problem would
become very complex and is not attempted here. It seems reasonable to conjecture
that if if oe 0 2
will suffice to solve the linear equations for
the correction w (k) using the Cholesky factorization of
I). Methods
for the solution of the linear systems are considered in more detail in Section 4.
4 Solving the linear systems.
In the RQI method formulated in the previous section the main work consists
of solving in each step two linear systems of the form
Here oe is an approximation to oe n+1 and varies from step to step. Provided that
the system (4.1) is symmetric and positive definite.
4.1 Direct linear solvers.
then the system (4.1) can be solved by computing the (sparse)
Cholesky factorization of the matrix A T A \Gamma oe 2 I. Note that A T A only has to
be formed once and the symbolic phase of the factorization does not have to be
repeated. However, it is a big disadvantage that a new numerical factorization
has to be computed at each step of the RQI algorithm.
For greater accuracy and stability in solving LS problems it is often preferred
to use a QR factorization instead of a Cholesky factorization. However, since
in the TLS normal equations the term oe 2 I is subtracted from A T A, this is not
straightforward. The Cholesky factor of the matrix A T A \Gamma oe 2 I can be obtained
from the QR factorization of the matrix
A
ioeI
, where i is the imaginary unit.
This is a downdating problem for the QR factorization and can be performed
using stabilized hyperbolic rotations, see [2, pp. 143-144], or hyperbolic Householder
transformations, see [22]. However, in the sparse case this is not an
attractive alternative, since it would require nontrivial modifications of existing
software for sparse QR factorization.
4.2 Iterated deregularization.
To solve the TLS normal equations using only a single factorization of A T A
we can adapt an iterated regularization scheme due to Riley and analyzed by
Golub [9]. In this scheme, we solve the TLS normal equations by the iteration
A T Affi
If lim k!1 x b. This iteration will converge with
linear rate equal to ae
n provided that ae ! 1. This iteration may be
implemented very efficiently if the QR decomposition of A is available. We do
not pursue this method further, since it has no advantage over the preconditioned
conjugate gradient method developed in [3].
4.3 A preconditioned conjugate gradient algorithm.
Performing the change of variables is a given nonsingular
matrix, and multiplying from the left with S \GammaT the system (4.1) becomes
This system is symmetric positive definite provided that oe ! oe 0 n , and hence
the conjugate gradient method can be applied. We can use for S the same
preconditioners as have been developed for the LS problem; for a survey see [2,
Ch. 7].
In the following we consider a special choice of preconditioner, the complete
Cholesky factor R of A T A (or R from a QR decomposition of A). Unless A is
huge this is often a feasible choice, since efficient software for sparse Cholesky
and sparse QR factorization are readily available [2, Ch. 7]. Using AR
I, the preconditioned system (4.2) simplifies to
(Note that although A and A T have disappeared from this system of equations
matrix-vector multiplications with these matrices are used to compute the right-hand
side f !) In the inverse iteration step used in the initialization,
the solution obtained by two triangular solves.
The standard conjugate gradient method applied to the system (4.2) can be
formulated in terms of the original variables w. The resulting algorithm is a
slightly simplified version of the algorithm PCGTLS given in [3] and can be
Algorithm 4.1. PCGTLS
Preconditioned gradient method for solving using the
Cholesky factor R of A T A as preconditioner.
Initialize: w ks
.
For
ks (j+1) k 2fi
Denote the original and the preconditioned matrix by I and
e
respectively. Then a simple calculation shows that for
the condition number of the transformed system is reduced by a
factor of -(A),
!/
The spectrum of ~
C will be clustered close to 1. In particular in the limit when
the eigenvalues of e
C will lie in the interval
(Note the relation to the condition number - TLS !) Hence, unless oe 0
can expect this choice of preconditioner to work very well for solving the shifted
system (4.1).
The I is positive definite if oe ! oe 0 n . In this case
in PCGTLS, and the division in computing ff k can always be carried out. If
n then the system (4.2) is not positive definite and a division by zero
can occur. This can be avoided by including a test to ensure that
equivalently kp the CG iterations are considered to
have failed. The RQI step is then repeated with a new smaller value of oe 2
e.g.,
The accuracy of TLS solutions computed by Rayleigh Quotient Iteration will
basically depend on the accuracy residuals and the stability of the method used
to solve the linear systems (4.1). We note that the cg method CGLS1 for the
LS problem, which is related to PCGTLS, has been shown to have very good
numerical stability properties, see [5].
4.4 Termination criteria in PCGTLS.
The RQI iteration, using PCGTLS as an inner iteration for solving the linear
systems, is an inexact Newton method for solving a system of nonlinear equa-
tions. Such methods have been studied by Dembo, Eisenstat, and Steihaug [7],
who consider the problem of how to terminate the iterative solver so that the
rate of convergence of the outer Newton method is preserved.
Consider the iteration
where r k is the residual error. In [7] it is shown that maintaining a convergence
order of 1 requires that when k !1, the residuals satisfy inequalities
is a forcing sequence.
In practice the above asymptotic result turns out to be of little practical use
in our context. Once the asymptotic cubic convergence is realized, the ultimate
accuracy possible in double precision already has been achieved. A more prac-
tical, ad hoc termination criterion for the PCGTLS iterations will be described
together with the numerical results reported below.
Remark. In the second linear system to be solved in RQI,
the right-hand side converges to x TLS . Hence it is tempting to use the value of
u obtained from the last RQI to initialize PCGTLS in the next step. However,
our experience is that this slows down the convergence compared to initializing
u to zero.
5 Numerical results.
5.1 Accuracy and termination criteria.
Numerical tests were performed in Matlab on a SUN SPARC station 10 using
double precision with unit roundoff . For the initial testing we used
contrived test problems [A; similar to those in [5] and generated
in the following way. 1 Let
e
where Y; Z are random orthogonal matrices and
Further, let
Ax:
This ensures that the norm of the solution does not reflect the illconditioning of
A. We then add random perturbations
Note that since oe there is a perturbation E to A with
which makes A rank deficient. Therefore it is not realistic to consider perturbations
with To test the termination criteria for the inner iterations
iterations
log0||x

Figure

5.1: Errors problem PS(30,15), with
systems solved by PCGTLS with iterations.
iterations
log0||x

Figure

5.2: Errors
systems solved by PCGTLS with
we used problem P (30; 15), oe 0
The linear systems arising in RQI were solved using PCGTLS with the Cholesky
factor of A T A as preconditioning. The criterion (4.6) shows that the linear systems
should be solved more and more accurately as the RQI method converges.
The rate of convergence depends on the ratio oe n+1 =oe 0 n , see (4.4), and is usually
very rapid. We have used a very simple strategy where in the kth step of RQI
These test problems are neither large nor sparse!
iterations are performed, where - 0 is a parameter to be
chosen.
In

Figure

5.1 we show results for 2. The plots for
are almost indistinguishable, whereas delay in convergence.
Indeed, for this problem taking iterations in PCGTLS suffices to give
the same result as using an exact (direct) solver. Since no refactorizations are
performed the object should be to minimize the total number of PCGTLS iter-
ations. Based on these considerations and the test results we recommend taking
should work well for problems where the ratio oe n+1 =oe 0 n
is smaller.
Rarely more than 2-3 RQI iterations will be needed. In Figure 5.2 we show
results for problem PS(30,15), and different error levels
Here respectively, were needed to achieve an accuracy
of about 10 \Gamma11 in x TLS . Since oe 0 this is equal to the
best limiting accuracy that can be expected. Note also that the error in oe n+1
converges to machine precision, usually in one less iteration, which supports the
use of the criterion (3.17) to terminate RQI.
5.2 Improvement from inverse iteration.
We now show the improvement resulting from including an initial step of
inverse iteration. In Figure 5.3 we show results for the problem considered
above. For the first two error levels only one RQI iteration now suffices. For the
highest error level oe n+1 converges in two iterations and x TLS in three.
-22
iterations
log0||x

Figure

5.3: Errors
One step of inverse iteration. Linear systems solved by PCGTLS with k +1
iteration.
We now consider the second test problem in [14], which is defined
where A 2 R n\Thetan\Gamma1 . Here e is a vector with entries generated randomly from a
normal distribution with mean 0:0 and variance 1:0, and scaled so that jjejj
. For 0:01 the condition
numbers in (2.5)-(2.7) are
respectively. This problem has features similar to those of the small illcondi-
tioned example discussed previously in Section 2.2, although here the norm of
the solution x LS is large.
-22
iterations
log0||x

Figure

5.4: Second test problem with 0:001. RQI without/with one step of inverse
iteration
Applying the RQI algorithm we obtained the results shown in Figure 5.4. The
initial approximation ae(x LS ) is here far outside the interval [oe n+1 ; oe 0 n ). Thus
the matrix A T A \Gamma oe 2 I is initially not positive definite and we cannot guarantee
the existence of the Cholesky factor. However, the Algorithm PCGTLS still does
not break down, and as shown in Figure 5.4 the limiting accuracy is obtained
after five RQI iterations. This surprisingly good performance of RQI can be
explained by the fact that even though x LS does not approximate x TLS well,
the angle between them is small; the cosine equals 0:98453.
Performing one step of inverse iteration before applying the RQI algorithm
gives much improved convergence. The one initial step of inverse iteration here
suffices to give an initial approximation in the interval [oe n+1 ; oe 0 n ). This can
be compared with 12-23 steps of bisection needed to achieve such a starting
approximation, see [14]! Three RQI iterations now give the solution x TLS with
an error close to the limiting accuracy, see Fig. 5.4.
We note that in both cases we obtained oe n+1 to full machine precision. Also,
the relative error norm of in the TLS solution was consistent with the condition
5.3 A problem in signal restoration.
The Toeplitz matrix used in this example comes from an application in signal
restoration, see [14, Example 3]. Specifically, an n \Theta (n \Gamma 2!) convolution matrix
T is constructed to have entries in the first column given by
exp
and zero otherwise. Entries in the first row given by t
zero otherwise, where 8. A Toeplitz matrix T and right-hand
side vector g is then constructed as
e, where E is a
random Toeplitz matrix with the same structure as T , and e is a random vector.
The entries in E and e are generated randomly from a normal distribution with
mean 0.0 and variance 1.0, and scaled so that
In [14] problems with convergence were reported. However, these are due to the
choice of right-hand side -
1 , which was taken to be a vector of all ones. For the
unperturbed problem this vector is orthogonal to the space spanned by
the left singular vector corresponding to the smallest singular value. Therefore
the magnitude of the component in this direction of the initial vector x LS will be
very small, of the order fl. Also, although A is quite well conditioned the least
squares residual is large. The TLS problem is therefore close to a nongeneric
problem and thus very illconditoned.
Because of the extreme illconditioning for this right-hand side, the behavior
of any solution method becomes very sensitive to the particular random perturbation
added. We have therefore instead chosen a right-hand side - g 2 given
by - m. For this the TLS problem is much better
conditioned, see Table 5.1. Convergence is now obtained in just two iterations,
see

Figure

5.5.
Table

5.1: Condition numbers for test problem 3 for right-hand sides - g i ,
-2iterations
log0||x

Figure

5.5: Third test problem; RQI with one step of inverse iteration,
6

Summary

.
We have developed an algorithm for solving large scale TLS problems based
on Rayleigh quotient iteration for computing the right singular vector of
defining the solution. The main work in this method consists of solving a sequence
of linear systems with matrix A T A \Gamma oe 2 I, where oe is the current approximation
to the smallest singular value of oe n+1 of (A; b). For large and sparse
TLS problems these linear systems can be solved by a preconditioned conjugate
gradient method. An efficient preconditioner is given by a (possibly incomplete)
Cholesky factorization of A T A or QR factorization of A.
Termination criteria for the inner and outer iterations have been given. We
conjecture that the described method almost always computes the TLS solution
with an accuracy compatible with a backward stable method. Although
a detailed error analysis is not given this conjecture is supported by numerical
results.
Methods for solving the TLS problem are by necessity more complex than
those for the (linear) LS problem. Our algorithm contains several ad hoc choices.
On the limited set of test problems we have tried it has only failed for almost
singular problems, for which the total least squares model is not relevant and
should not be used.
In our method the perturbation E is a rank one matrix which in general is
dense. Sometimes it is desired to find a perturbation E that preserves the
sparsity structure of A. A Newton method for this more difficult problem has
been developed by Rosen, Park, and Glick [21]. However, the complexity of this
algorithm limits it to fairly small sized problems. Recently a method, which has
the potential to be applied to large sparse problems has been given by Yalamov
and Yun Yuan [26]. Their algorithm only converges with linear rate, which may
suffice to obtain a low accuracy solution.



--R








Improving the accuracy of computed singular values.
Numerical methods for solving least squares problems
An analysis of the total least squares problem

A survey of condition number estimation for triangular matrices
Accuracy and Stability of Numerical Algorithms
A total least squares method for Toeplitz systems of equations
The minimum eigenvalue of a symmetric positive-definite Toeplitz matrix and rational Hermitian interpolation
Sparse QR factorization in MATLAB

The Symmetric Eigenvalue Problem
On the convergence of a practical QR algorithm

Total least norm formulation and solution for structured problems
stability and pivoting
Criteria for combining inverse and Rayleigh quotient iteration
Iterative algorithms for computing the singular subspace of a matrix associated with its smallest singular values
The Total Least Squares Problem: Computational Aspects and Analysis
A successive least squares method for structured total least squares.
Iterative Methods for Least Squares and Total Least Squares Problems
--TR

--CTR
Computing smallest singular triplets with implicitly restarted Lanczos bidiagonalization, Applied Numerical Mathematics, v.49 n.1, p.39-61, April 2004
rayleigh quotient iteration;singular values;conjugate gradient method;total least squares