Abstract
The constrained least-squares n × n-matrix problem where the feasibility set is the subspace of the Toeplitz matrices is analyzed. The general, the upper and lower triangular cases are solved by making use of the singular value decomposition. For the symmetric case, an algorithm based on the alternate projection method is proposed. The implementation does not require the calculation of the eigenvalue of a matrix and still guarantees convergence. Encouraging preliminary results are discussed.
constrained least squares; Toeplitz matrix; singular value decomposition; alternate projection method
Finding the closest Toeplitz matrix* * This work has been supported by UNS grant 24/L031.
Maria Gabriela Eberle; Maria Cristina Maciel
Department of Mathematics, Southern National University Av. Alem 1253, 8000 Bahía Blanca, Argentina
E-mail: geberle@criba.edu.ar, immaciel@criba.edu.ar
ABSTRACT
The constrained least-squares n × n-matrix problem where the feasibility set is the subspace of the Toeplitz matrices is analyzed. The general, the upper and lower triangular cases are solved by making use of the singular value decomposition. For the symmetric case, an algorithm based on the alternate projection method is proposed. The implementation does not require the calculation of the eigenvalue of a matrix and still guarantees convergence. Encouraging preliminary results are discussed.
Mathematical subject classification: 65F20, 65F30, 65H10.
Key words: constrained least squares, Toeplitz matrix, singular value decomposition, alternate projection method.
1 Introduction
In this work the constrained least-squares Toeplitz matrix is analyzed. Let us begin considering the general constrained least squares matrix problem
where A, B Î mxn, with m > n and Ì nxn has a particular pattern. Throughout this paper A
denotes the Frobenius norm of A Î mxn defined as
In the last two decades an extensive research activity has been done on this kind of problems which arise in many areas. For instance, in statistics, the problem of finding the nearest symmetric positive definite patterned matrix to a sample covariance matrix is a classical example, [11]. The symmetric positive semidefinite case has analyzed by Higham [9].
In the investigation of elastic structures, the determination of the strain matrix which relates forces fi and displacement di according to Xfi = di is stated as (1). The orthogonal case has been studied by Schonemann [13] and the Higham [8] and the symmetric case by Higham [10]. These problems are known in the literature as the orthogonal and symmetric Procrustes problems respectively. When is the subspace of persymmetric matrices (that is, it is symmetric about the NE-SW diagonal) the resulting problem has been studied by Eberle [4, 5].
Now in this work we are interested in solving the constrained least squares Toeplitz matrix problem
where Ì nxn is the subspace of Toeplitz matrices.
Let us recall that a matrix T = (tij) Î nxn is said to be a Toeplitz matrix if it is constant along its diagonals, i.e there exist scalars , rn+1, ¼, r0,¼, rn1, such that tij = rji, for all i,j:
The rest of the paper is organized as follows. In section 2, the general Toeplitz problem is studied. The special cases of triangular and symmetric Toeplitz problems are analyzed in section 3 and 4 respectively. The numerical results are presented in section 5 and our conclusions and final remarks in section 6.
2 The general problem
In this section the feasibility set, , of problem (2) is characterized.
Given the matrices A, B Î mxn, the subspace can be represented by
where ap Î and Gp Î nxn are matrices defined as follows
We must observe that the matrices Gp are characterized by the following fact: for each entry ij, there exists an unique index p, (n1) < p < (n1), such that (Gp)ij = 1. Clearly they form a basis for the subspace .
2.1 The transformed problem
We will transform the problem (2) in a simpler one. To do this transformation let A be the matrix have the following singular value decomposition
where P Î mxm and Q Î nxn, are orthogonal matrices and S = diag(s1,¼, sn), s1> ¼ > sn > 0.
Since the Frobenius norm remains invariant under orthogonal transformations it is obtained
with
then the problem 2 is equivalent to
where ¢ is the following subspace:
2.2 Characterization of the solution on ¢
We introduce the notation P(A) meaning the projection of the matrix A on the set . The following theorem characterizes the projection on the subspace of matrices ¢.
Theorem 2.1. If C1Î nxn, then the unique solution of the problem (4) is
for all (n1) < l < (n1).
Proof. The objective function f : 2n1
is given by
Since f is twice continuously differentiable, we compute (a) for all p such that (n1) < p < (n1)
The factor in the right hand side of (5) can be written as
and (5) becomes in
Having in mind that
is the inner product between the ith row of (SQT) and the jth column of alGl, we obtain
Similarly, the element (SQTGp)ij, is the inner product between the ith row of SQT and the column j of Gp. Then:
(SQTGp)ij = siQ(jp)i .
Because of the sparse structure of Gp the entries (Gp)ij = 1 if j > p.
If p > 0, then (Gp)ij = 1 for all j such that p + 1 < j < n.
If p < 0, then (Gp)ij = 1 for all j such that 1 < j < (n + p).
Therefore (6) can be written as
Defining
with
vT = (Q(jk)1Q(jp)1, Q(jk)2Q(jp)2,¼,Q(jk)nQ(jp)n),
q = ,
we have that the sum of all the components of v is the inner product between two columns of the orthogonal matrix Q and then it is zero when k ¹ p. Therefore for k ¹ p we have:
then vTq = 0 for all k ¹ p and
Then (6) becomes
for all p = 1,¼, n. Therefore, from the first order necessary condition
we obtain
We observe that the denominator of (7) is non zero, because if it were zero for any i, j, it would be
Now, we need to verify that is a minimizer of f. We compute
and
It implies that the Hessian matrix Ñ2f() is positive define and therefore is the minimizer of f.
3 The triangular Toeplitz problem
In this section we consider the case when the feasibility is the subspace of triangular Toeplitz matrices. We begin considering the upper triangular case.
3.1 The upper triangular Toeplitz problem
Let the minimization problem be
where A, B Î mxn and is the subspace of upper triangular Toeplitz matrices. If X Î , then
The subspace , can be written as:
where apÎ and the matrices Gp Î nxn are defined for all 1 < i, j, p < n as follow
We remark that the matrices Gp are characterized by the property that for each entry ij, there exists an unique p, 1 < p < n such that (Gp)ij = 1.
The set , is a basis for the subspace of the upper triangular Toeplitz matrices of nxn. Similarly to the general Toeplitz case, the singular value decomposition is used in order to obtain a solution of an equivalent problem
with
P Î mxm y Q Î nxn, are orthogonal matrices, S = diag(s1,¼, sn) with s1> ¼ > sn> 0, and the subspace
3.1.1 Characterization of the solution on ¢
The following theorem, similar to Theorem 2.1 for the dense case characterizes the projection on ¢.
Theorem 3.1. Let be C1Î nxn, then the unique solution of the problem
is given by
for all 1< l < n.
Proof. Similar to the proof of theorem 2.1 having in mind that the first derivative is
The element (SQT((alGl))ij is the inner product between the ith row of (SQT) and the jth column of (alGl), then it results
Analogously, the element ( SQTGp)ij, is the inner between the row i of SQT and the column j of Gp, then
3.2 The lower triangular Toeplitz problem
Let us consider the problem
where A, B Î mxn and is the subspace of the lower triangular Toeplitz matrices. If X Î , then
The subspace , can be expressed as
where bpÎ and the matrices Gp Î nxn are defined as
for all 1 < i, j, p < n. The set , is a basis for the subspace . Again, applying the singular value decomposition of A the problem can be transformed in
with
where P Î mxm y Q Î nxn, are the orthogonal matrices, S = diag(s1,¼, sn), s1> ¼ > sn > 0, and the subspace
As in the general and upper triangular Toeplitz cases the projection is characterized on ¢.
Theorem 3.2. If C1Î nxn, then the unique solution of the problem
is given by
for all 1 < l < n.
Proof. Similar to the proof of Teorema 2.1.
4 The symmetric Toeplitz problem
This section is concerned with the problem
where A, B Î mxn , is the subspace of the Toeplitz matrices and is the subspace of the symmetric matrices of nxn. Clearly, the solution is a matrix which is symmetric not only with respect to the main diagonal but also to the SW-NE diagonal. Since the feasibility region is the intersection of two subspace and because we know how to solve the persymmetric Procrustes problem [5] and the general Toeplitz problem, see section 2, we develop an algorithm based on the alternating projection method [3, 1, 6, 7]. Therefore the solution of (16) is obtained projecting alternating on the subspaces y .
Given the problem
according with Escalante and Raydán [7] it can be transformed, via the singular value decomposition of A in the following equivalent problem
where
and the feasibility region
¢ = {Z Î nxn : Z = SY, Y = YT}.
Hence if the solution of the last problem is
where is obtained by the methods proposed by Higham for the symmetric case [10], the solution of the first problem is
On the other hand, we know how to solve the general Toeplitz problem, see section 2. Combining both results, to solve the problem 16 is equivalent to solve
where
Remark. The set ¢¢, comes from the fact that the symmetric and Toeplitz problems have the same objective function.
The alternating projection algorithm for the symmetric Toeplitz case (19) is as follows:
Algorithm 4.1. Given A, B Î mxn. X0 = C1Î nxn, where C1 is obtained via the singular value decomposition of A, for i = 0, 1, 2, ¼, until convergence repeat
The algorithm terminates when two consecutive projection on the same subspace are close enough. The convergence is guaranteed by the general theory established by von Neumann [12] for the alternating projection algorithm.
5 Numerical experience
In this section we show some examples. The algorithms have been implemented in Fortran, by using FORTRAN POWER STATION (1993) in an environment PC with a processor Pentium(R2) Intelmmx(TM) Technology, 31.0 MB de RAM. and 32 bits of virtual memory. The LINPACK library dsvdc subroutine have been used to obtain the singular value decomposition.
Example 5.1. The general problem (2) has been solved with the following matrices
The solution is
We project on the sets y .
a) The projection on is
b) The projection on is
Example 5.2. We illustrate the behavior of the algorithm 4.1 for the symmetric case carrying out the following example
The solution is
The algorithm terminates when two consecutive projections on the subspace of the Toeplitz matrices is less or equal to a fixed tolerance. We set = 10-6 for this experiment.
6 Concluding remarks
We have studied the Procrustes Toeplitz problem, analyzing the general, the triangular and the symmetric cases. In the last case an algorithm based on the alternating projection method is proposed. We illustrate our results with some examples.
The Toeplitz systems Tx = b where T Î nxn is a Toeplitz matrix, are present in many disciplines (signal and image processing, times series analyzes, queuing networks, partial differential equations problems and control problems). When these systems are solved via preconditioned conjugate gradient methods, the search of adequate preconditioners is the main issue. Among the possible preconditioners for a Toeplitz matrix, Chan proposes the preconditioner c(T) defined to be the minimizer of C T
. Therefore a possible direction of future work is the application of the results obtained to answer questions related to Toeplitz systems. The reader interested in this subject can review the survey of Chan and Ng [2] and the references therein.7 Acknowledgments
The authors wish to thank two anonymous referees for helpful suggestions.
#537/01.
- [1] J.P. Boyle and R.L. Dykstra, A method for finding projections onto the intersections of convex sets in Hilbert spaces, Lecture Notes in Statistics, 37 (1986), 28-47.
- [2] R.H. Chan and M.K. Ng, Conjugate gradient methods for Toeplitz systems, SIAM Review, 38(3) (1996), 427-482.
- [3] R.L. Dykstra, An algorithm for restricted least squares regression, Journal of the Amer. Stat. Ass., 78(384) (1983), 837-842.
- [4] M.G. Eberle, Métodos de proyecciones alternas en problemas de optimización persimétricos, Master's thesis, Dept. de Matemática, Universidad Nacional del Sur, Bahía Blanca, Argentina, (1999).
- [5] M.G. Eberle and M.C. Maciel, The persymmetric Procrusto problem, Bahía Blanca, Argentina, (2001). In preparation.
- [6] R. Escalante and M. Raydan, Dykstra's algorithm for constrained least-squares matrix problem, Num. Lin. Alg, with Appl., 36 (1996), 459-471.
- [7] R. Escalante and M. Raydan, Dykstra's algorithm for constrained least-squares rectangular matrix problem, Computers and Mathematics with Applications, 35(6) (1998), 73-79.
- [8] N.J. Higham, Newton's methods for the matrix square root, Mathematics of Computation, 46(174) (1986), 537-549.
- [9] N.J. Higham, Computing a nearest symmetric positive semidefinite matrix, Linear Algebra and its Applications, 103 (1988), 103-118.
- [10] N.J. Higham, The symmetric Procrustes problem, BIT, 28 (1988), 133-143.
- [11] H. Hu and I. Olkin, A numerical procedure for finding the positive definite matrix closest to a patterned matrix, Statistics and Probability Letters, 12 (1991), 511-515.
- [12] J. Von Neumann, Functional operator, vol. ii: The geometry of orthogonal spaces, In Annals of Math. Studies. Princeton University Press, Princeton, (1950).
- [13] P.H. Schonemann, A generalized solution of the orthogonal Procrustes problem, Psychometrica, 31 (1966), 1-10.
Publication Dates
-
Publication in this collection
19 July 2004 -
Date of issue
2003