Fast Algorithms for Structured Least Squares and Total Least Squares Problems.

We consider the problem of solving least squares problems involving a matrix M of small displacement rank with respect to two matrices Z 1 and Z 2. We develop formulas for the generators of the matrix M (H) M in terms of the generators of M and show that the Cholesky factorization of the matrix M (H) M can be computed quickly if Z 1 is close to unitary and Z 2 is triangular and nilpotent. These conditions are satisfied for several classes of matrices, including Toeplitz, block Toeplitz, Hankel, and block Hankel, and for matrices whose blocks have such structure. Fast Cholesky factorization enables fast solution of least squares problems, total least squares problems, and regularized total least squares problems involving these classes of matrices.

1. Introduction

In [6], Mastronardi, Lemmerling, and Van Huffel present an algorithm for solving fast structured total least squares problems of the form subject to the constraints with ∈ C× a given Toeplitz matrix and ∈ C×1 a given vector. They include one additional constraint: is a Toeplitz matrix. They produced a fast algorithm for solving this structured total least squares problem (STLS) and showed that the solution was a better estimator than the solution to the total least squares problem without the Toeplitz constraint.

In this paper, we work through the details of a small generalization of this algorithm. We consider the same problem (1), but under the constraint that and have small displacement rank relative to some matrices 1 and 2. Choosing these two matrices to be shift-down matrices and the rank to be two gives the Toeplitz constraint considered by [6], but we will be interested in other cases as well. See [3] for a discussion of displacement structure.

We also consider fast solution of the problem under the additional constraint that the norm of the solution vector is specified, a problem posed in Pruessner and O’Leary [10]. This corresponds to a Tikhonov regularization of our structured total least squares problem and results in a fast solution algorithm for the problem considered in [8, 9, 10, 7].

The core of the algorithm in [6], based on a more general algorithm of [11], relies on two results: the representation of the generators for the matrix that appears in the normal equations when is Toeplitz, and then a fast factorization of a matrix derived from these generators. So we begin in Sec. 2 with a construction of the generators for when is any matrix of small displacement rank and 1 is close to unitary. In Sec. 3 we show that it is inexpensive to form a Cholesky factorization of whenever 2 is triangular and nilpotent. Section 4 concerns the application of this algorithm, and we conclude in Sec. 5.

The results in this work are derived from [4, 5].

2. Generators of M

One way to solve a least squares problem involving a matrix of full column rank is to use the normal equation formulation; we factor = where is a square lower triangular matrix and then solve the linear system

Thus, we can solve this problem fast if we can compute a Cholesky factorization of fast. To do this, we will first derive a generator for the matrix when ∈ C× (m ≥ n) has low displacement rank.

Suppose that has low displacement rank relative to the matrices 1 ∈ C× and 2 ∈ C×, which means that if we define then rank (), which we denote by ρ1, is small relative to n.

Suppose is a unitary matrix , where has rank ρ2, also assumed to be small. For example, if is Toeplitz, let 1 be the shift-down matrix with ones on its subdiagonal and zeros elsewhere, and then is the matrix with a one in the last position of row 1.

Then also has low displacement rank relative to 2, as we can see from the identity

Theorem 1. If the rank of is ρ1 and if the unitary matrix is equal to 1 + where has rank ρ2, then has rank at most 2(ρ1 + ρ2).

Proof: The equation in the statement of the theorem is a regrouping of the terms in the previous equation. The rank of is at most the rank of plus the rank of , so the rank of the sum in that equation is at most 2(ρ1 + ρ2). []

An alternate proof of this theorem can be derived by noting that displacement rank is preserved by taking a Schur complement [3, Lemma 1.5.2], and is the Schur complement of which has displacement rank at most 2(ρ1 + ρ2) relative to the matrix

The bound 2(ρ1 + ρ2) can be an overestimate, since the and terms may interact. For example, for a Toeplitz matrix with rank computed relative to the shift-down matrices, the bound yields 6 while the actual rank is 4. If we use circular shifts, though, both the bound and the actual rank are 4.

Using the identity

3. Determining a Factorization of M From the Generators

Theorem 1 tells us how to determine ρ = 2 (ρ1 + ρ2) vectors so that where s equals plus or minus 1. When 1 and 2 are shift-down matrices, it has been shown [6, 1, 2] that this implies that where = diag(), = diag(s) (n × n), and () is the lower triangular Toeplitz matrix with first column equal to . We now generalize this result somewhat.

3.1 Triangular Factors of Structured Matrices

Theorem 2.

If 1 is nilpotent, then and only if where

Proof: Suppose . Observe that and so, since , we conclude that

To prove the converse, suppose . Then, since we conclude that if , then

Now since 1 is nilpotent, for some p ≤ n. Therefore, , and working backward in powers of 1, we see that , so . []

The following corollary can be proved by finite induction.

Corollary 1: If 1 is nilpotent, then if and only if

3.2 Triangular Factors of M

In order to solve our least squares problem, we wish to determine a Cholesky factorization so we need to reduce the matrix to upper triangular form, where the vectors are defined in (3).

If 1 and 2 are shift-down matrices, then [6] shows how to do this reduction fast. Using Corollary 1, we will see that this can be done fast whenever 2 is triangular and nilpotent. We present the algorithm for the case 2 lower triangular; the upper triangular case is analogous but works with the columns in reverse order.

The algorithm proceeds by columns, putting zeros below the main diagonal. Note that

Suppose we determine a rotation between the first row and row n + 1, which contains , to zero the first element of . The same rotation between and also zeroes the first element of since is upper triangular. Therefore, by introducing one zero into our matrix, we have implicitly introduced n − 1 more, so we can put zeroes below the main diagonal in column 1 by using only ρ − 1 rotations, independent of the size of n.

We then use the resulting second row, equal to the first row postmultiplied by , to zero the second element of row n + 1. Again this implicitly introduces additional zeros, n − 2 of them, and we complete the operations on column 2 by using ρ − 1 rotations.

If we repeat this for each column, we accomplish our reduction. Let be the ρ × n matrix whose rows are . We can thus reduce to upper triangular form just by operating on the matrix .

We design our algorithm to use Givens rotations as often as possible, minimizing the number of hyperbolic rotations in order to preserve stability. A Givens rotation can be used between row i and row j whenever s and s (see equation (3)) have the same sign; if they have different signs, then we must use a hyperbolic rotation. We’ll assume that we have ordered the generators so that the first rows of have s = 1 and the remaining ones have s = − 1.

The cost of this reduction is at most O(ρn2), ignoring sparsity, plus the cost of the multiplications by 2. Sometimes sparsity can reduce this cost significantly [7]. Without exploiting the structure of the cost would be O(ρn3). Once the factors are computed, they can then be used to solve (2).

4. Some Applications

Our fast algorithm enables us to solve least squares, total least squares, and regularized total least squares problems involving matrices for which 1 is close to unitary and 2 is triangular and nilpotent. This includes several important classes of structured matrices.

4.1 Solving Least Squares Problems

Consider first the least squares problem where ∈ C × has full column rank. We can use our algorithm if

is Toeplitz. Then 1 and 2 are the shift-down matrices of appropriate size and the displacement rank is 2.

is block Toeplitz: where has dimension m / γ × n / δ. Then choosing 1 ∈ C× as where each block has dimension m/γ×m/γ, and choosing 2∈Cm×n similarly but with blocks of dimension n / δ × n / δ gives a displacement rank of m / γ + n / δ.

has Toeplitz blocks: where each block is Toeplitz. Then choosing 1 ∈ C× as a block-diagonal matrix consisting of shift-down matrices of dimensions matching the column dimension of the diagonal blocks of and choosing 2 ∈ C× similarly but with blocks matching the row dimension gives a displacement rank of 2(γ + δ).

Hankel, block Hankel, or having Hankel blocks. These cases are analogous to those considered above.

4.2 Solving Structured Total Least Squares Problems

In some cases the matrix of the problem can be estimated only with error, and we need to determine not just the parameters of the least squares fit but also the corrections to the operator. This problem can be formulated as subject to the constraints with and having the same structure.

Suppose that the matrix can be specified by p parameters α1,…, α. For example, if E is a Toeplitz matrix, then and p = m + n − 1. We rewrite our problem as where

Following [6], we have replaced the term by , equivalent except for scaling of the entries .

We define the matrix ∈ C× by the equation

For example, if is Toeplitz, then p = m + n − 1 and

Following [11], we form a quadratic approximation to (5) by using linear approximations + and + , resulting in so that

If we minimize this with respect to and , then we can form a new approximation to the solution of (5) and then repeat the procedure until convergence. As noted by [11], this is a Gauss-Newton algorithm applied to (5). Therefore, the main computational task is to solve linear least squares problems of the form where

If is in one of the classes considered in Sec. 4.1, then the matrix has low displacement rank and we can solve the problem fast.

4.3 Solving Regularized Total Least Squares Problems

In many deblurring problems and other discretized problems involving integral equations of the first kind, the matrix is so ill-conditioned that noise in the observations is magnified in solving the STLS problem and a meaningful solution cannot be obtained.

In this case it is necessary to add a regularization constraint to the problem. One common regularization constraint is to restrict the size of the solution, or some linear transformation of the solution: where u is a given scalar and is commonly chosen to be the identity matrix or a difference operator. If has low displacement rank relative to 2 and an appropriate , our algorithm can be easily modified to incorporate regularization. In this case, our problem (5) can be reformulated as where = ( + ) – and λ, the regularization parameter, is the Lagrange multiplier for the new constraint. Using a derivation similar to that above, the linearization of (7) results in the following problem to be solved at each step of the iteration:

Thus, our new matrix is the matrix from the previous section augmented by the extra rows [0, λ], and the only change necessary in the algorithm is to find the generators of this matrix rather than the old one.

The displacement structure of this matrix is greatly simplified if is upper triangular and 1 and 2 are shift-down matrices. As noted before, is zero except for a one in the last position of the first row, and thus is zero except for a λ in the last position of the first row. Therefore, , so, applying Theorem 1, we have the following result.

Theorem 3. If is upper triangular and 1 and 2 are shift-down matrices, then and has rank 2ρ1, where ρ1 is the rank of . More generally, (8) holds whenever .

5. Conclusions

We have derived the generators for when is any matrix of small displacement rank relative to 1 and 2. We have shown that it is inexpensive to form a Cholesky factorization of whenever 1 is close to unitary and 2 is triangular and nilpotent, and we have generalized this algorithm when a regularization constraint is to be applied.

Algorithm Reduce ()3

For j = 1, …, n,

For i = 2, …, ρ^,

If hij is nonzero, then

zero it by a Givens rotation between row 1 and row i;

end for

For i=ρ^+2,…,ρ,

If hij is nonzero, then

zero it by a Givens rotation between row ρ^+1 and row i;

end for

If h_ρ^+1,j is nonzero, then

zero it by a stabilized hyperbolic rotation between row 1 and row ρ^+1;

Then the jth row of LH is h1H, the first row of the current H matrix.

Replace the first row of H by h1HZ2H to form the pivot row for the next value of j.

end for