Properties of the Toeplitz Covariance Matrix Numerical Likelihood Ratio Maximization: Convex-Like Optimization Behavior

Introduction

The Maximum Likelihood (ML) estimation of the Toeplitz covariance matrix, given a set of T independent identically distributed (i.i.d.) complex Gaussian N-variate vectors, is one of the classical signal processing problems that has been under intensive investigation since at least the ‘80s, and yet it continues to attract attention nowadays [3-19]. One of the main problems is that the nature of this optimization problem remains unclear. The type of this optimization problem was not specified in the seminal paper by J.P. Burg, D.G. Luenberger, and D.L. Wenger in 1982 or in the papers by D.B. Fuhrmann, but, written in 2024, the problem was described as non-convex [7,20]. Therefore, the nature of this important optimization problem remained unclear despite a large number of efficient solutions proposed over the decades [3-17]. Correspondingly, the properties of the globally optimal maximum likelihood Toeplitz covariance matrix estimates have not been explored because of the expectation of non-convexity, with potentially many solutions and the inability to specify the globally optimal ones.

In our recent papers, we addressed this problem by modifying the maximized likelihood function into a likelihood ratio with a probability density function (pdf) for the true covariance matrix, which does not depend on the true matrix and is specified by a priori known parameters [1,2]. The application of the standard MATLAB fmincon optimization routine revealed that while the overwhelming majority of solutions were non-positive definite and therefore easily rejected, a small number of positive definite solutions possessed a very important property. Specifically,irrespective of the initial solution of this iterative optimization routine that ended up in a positive definite (p.d.) solution, this solution was the same. Moreover, the maximized likelihood ratio exceeded that of the true Toeplitz covariance matrix. In fact, all acceptable (i.e., positive definite) solutions of this optimization problem possessed these most important properties of the globally optimum solution.

In this paper, we continue to investigate this phenomenon. Specifically, by applying the Carathéodory description of the optimized Toeplitz matrix, we avoided non-positive definite solutions but did not achieve the LR maximum obtained by direct covariance lag optimization. Therefore, after fmincon convergence on the Carathéodory parameters, we switched to direct optimization of the covariance lags. The overwhelming majority of the solutions possessed the properties of the LRoptimum p.d. Toeplitz matrix. Yet, the "distance" between this solution with its improved LR and true covariance matrix was proven to be greater than the p.d. solutions in Carathéodory parameterization. This "distance" was measured by the likelihood ratio with the sample covariance matrix replaced by the true Toeplitz covariance matrix .

Finally, we demonstrate that the application of this Toeplitz covariance matrix estimation to the problem of direction-of-arrival (DOA) estimation in a uniform linear antenna (ULA) array, allows it to significantly extend into the region of a smaller input SNR and/or smaller inter-source separation, the efficient resolution and DOA estimation of these sources, compared to MUSIC, applied to the sample covariance matrix .

By focusing on these new results, the paper concentrates on the properties of the solutions provided by the MATLAB fmincon routine and the gains in the classical DOA estimation problem these properties offer. Specifically, we demonstrate that the p.d. solutions, delivered by the MATLAB fmincon routine, possess the following properties:

• the optimized LR value exceeds the LR value produced by the true Toeplitz covariance matrix,

• the converged positive definite Toeplitz matrices are the same, irrespective of the initial solution used. These tested initial solutions included the very far from the true covariance matrix initiations with LR < 10−20, and the true covariance matrix .

These two properties belong to the truly globally optimal solution, but in our case, they describe only a subset of all possible solutions, including non-p.d. solutions. Since the non-p.d. solutions could be easily identified, the remaining p.d. solutions possess these important properties of the globally optimal ones.

The paper is focused on the two following goals:

• to describe a new technique for ML Toeplitz matrix estimation, which improves DOA estimation performance

• to investigate the properties of the quasi-global solution of the ML Toeplitz covariance matrix estimation. Our main contributions are summarized as follows.

i We demonstrate that while the MATLAB fmincon routine optimizes the covariance lags of the Toeplitz matrix, the vast majority of solutions are not positive definite.

ii All positive definite solutions provided by the MATLAB fmincon routine are practically the same, irrespective of the p.d. Toeplitz matrix used for initialization. The same solutions are obtained both with initializations that differ substantially from the true matrix (with LR ≤ 10^-20) and with initializations utilizing the true Toeplitz covariance matrix . The optimized LR values for these solutions exceed those produced by the true covariance Toeplitz matrix. The independence of the result from the initialization matrix, including the true covariance matrix, together with exceeding the LR value of the true covariance matrix, are the most important properties of the globally optimum solution. Yet, the numerical methods highlight these properties but do not provide analytical proof of global optimality.

iii For avoiding the non-p.d. solutions generated by the MATLAB fmincon routine, we proposed starting LR maximization using the Carathéodory p.d. Toeplitz matrix representation, followed by LR maximization with further covariance lag updates (i.e., in Vandermonde parameterization).

iv We demonstrated that this second stage of the LR optimization of the covariance lags increases the LR but decreases the proximity to the true covariance matrix that replaces the sample matrix in the likelihood ratio.

Correspondingly, in Sec 2, we provide the analytical description of the optimization problem and its possible parameterizations. In Sec 3, we demonstrate the properties of the optimized solutions for data with a Toeplitz matrix and describe the clutter returns in the HF OTHR. In Sec 4, we demonstrate that applying the proposed Toeplitz matrix estimation technique yields significant gains in source DOA estimation for sources impinging upon a uniform linear array. In Sec 5, we conclude our paper.

Maximum Likelihood Covariance Toeplitz Matrix Optimization Algorithm

Numerical Likelihood Ratio Maximization for Models with Different Eigenvalues of the Toeplitz Covariance Matrix

the covariance lag parameterization exceeded the LR from the Carathéodory parameterization alone by a few percent (up to 15%, Figure 1).

Figure 1: LR Gain Provided by the Final “Natural” Parameterization with Respect to the Initial “DoAs” Parameterization (1,000 trials)

Based on these numerical results, we may conclude that the proposed combination of parameterizations, starting with Carathéodory fmincon LR maximization and followed by the covariance lags parameterization, allowed us to mostly avoid solutions with negative eigenvalues, thus making an arbitrary p.d. Toeplitz matrix appropriate for fmincon initialization. By optimizing the covariance lags, we reach the LR maximum while mostly avoiding solutions with negative eigenvalues. Since the LR maximum is the same as the one achieved by fmincon when starting from the true covariance matrix , we may treat these solutions as globally optimum. Moreover, the number of these additional iterations for the covariance lag parameters is negligible.

In our first example (Figure 2), after 925 iterations with the Carathéodory parameters and LR = 1.708^.10−1, it required only 44 iterations using the Vandermonde covariance lags to reach the global LR maximum of LR = 1.853^.10−1.

Figure 2: LR Gain Provided by the Final “Natural” Parameterization with Respect to the Initial “DoAs” Parameterization (1,000 trials) In Figure 3, we show the percentage improvement in LR achieved by this additional optimization of covariance lags.

Figure 3: Equal LRs for Optimization, Started from the True Matrix (LR_Tn2opt) and a Random Matrix (LRopt) (1,000 trials)

One can see that in only one out of 1,000 conducted trials, the LR maximized over Carathéodory parameters exceeded the LR value from direct covariance lags optimization by no more than 0.01%. In all other trials, the direct optimization of covariance lags provided further LR improvement over the Carathéodory parameterization, with gains of less than 15% and an average of 6%. Therefore, starting the LR fmincon maximization from the Carathéodory matrix description, we avoid solutions with negative eigenvalues for a practically arbitrary p.d. Toeplitz matrix used to initialize the fmincon routine. The following LR maximization using the covariance lags description yields a globally optimal solution in all cases, with the achieved LR value equal to that initiated by the true covariance matrix .

Let us remind the reader that the comparative analysis of two different problem formulations in covariance lags and Carathéodory parameters was performed for the same sample covariance matrix . For this reason, in our next simulation series, we investigated 1,000 different sample matrices generated from the same data with the same true Toeplitz covariance matrix. The results of these simulations are illustrated in Figure 4 and Figure 5, where the dependence of the maximized LR on the number of iterations is provided.

Figure 4: fmincon LR Progression During the Initial Stage of Optimization in the “DoAs” Parameterization, Followed by the Optimization in “Natural” Parameterization

Figure 5: Another Example of the Two-Stage fmincon LR Maximization

The first part of the curve in Figure 4 corresponds to LR maximizations performed in Carathéodory parameters (925 iterations), followed by 44 iterations using covariance lag optimization by the fmincon routine. From ~150 iterations to the final 925 in the Carathéodory parameters, the LR value did not change much, finally reaching LR = 1.708 ^. 10−1. The next 44 iterations in covariance lags values, conducted by covariance lags trimming, reached the global maximum LR = 1.853^.10−1, which significantly exceeded the LR value produced by the true covariance matrix, LR = 1.494 ^. 10−1.

The second example, illustrated by the Figure 5, differs only in LR values, with the minimal gain provided by the second level LR maximization, reaching LR = 1.671 ^. 10−1 compared with the first-order LR optimization in Carathéodory parameters, reaching LR = 1.648 ^. 10−1. Note that the LR values generated by the true covariance matrix were equal to LR = 1.442 ^. 10−1, versus LR = 1.494^.10−1 in our first example.

The sample pdf of the LR gains (in percentage) achieved by the second-level LR maximization using the covariance lags optimization, compared with the LR achieved as a result of the first-level LR maximization in the Carathéodory parameters, is presented in Figure 6.

Figure 6: Distribution of the LR Gains Provided by the Second Stage LR Maximization in “Natural” Parameters

One can see that this gain can reach 10%, with an average of 3.2%. Only in very rare cases did the combined LR maximization slightly exceed the LR values obtained with the covariance lags parameterization alone.

These improvements do not exceed 0.005% and are therefore attributed to finite calculation accuracy. It seems acceptable that for quite complicated LR optimization algorithms, we may treat the results of our two-step LR maximization as globally optimal values, despite the negligible discrepancy caused by finite computational accuracy. In Figure 7, we provide the results of LR optimization, sorted in increasing order of the first step of LR maximization in Carathéodory parameters.

Figure 7: Maximized LR Values by “DoAs” Parameterization (blue line) and by the “Natural” Parameters Optimization (red line)

Then, for every LR value obtained from the first-step LR maximization over Carathéodory parameters, we provide the LR value achieved by the second-step LR maximization via direct trimming of the covariance lags. One can see that the initial LR maximization over Carathéodory parameters ranges from 0.12 to 0.28, while the second-stage LR maximization using the covariance lags adds, on average, 0.01 to the global LR maximum.

The conducted trials provide a detailed enough description of the statistical nature of the LR maximum using the proposed two-step numerical optimization.

C. Second Stage LR Maximization of the Covariance Lags Improves LR Values but Degrades the Proximity to the True Toeplitz Covariance Matrix

In statistical signal processing, it is well known that the more specific the a priori knowledge about the estimated set, the smaller the excess optimized likelihood ratio is and the closer the estimation result is to the true parameters. The comparison of the Hermitian and Toeplitz Hermitian covariance matrices' ML estimates is one such example, where the maximized LR value for the Hermitian matrix is always 1, whereas a priori knowledge of its Toeplitz structure leads to a significant degradation in the maximized LR value.

The discrepancy between the optimized LR values for the Carathéodory parameters and the covariance lag parameters should not arise, since we are only interested in positive definite matrices.

Yet, the demonstrated superior performance of the fmincon LR optimization for covariance lag parameters is beyond any doubt and is not associated with a more specific (more accurate) model of the estimated Toeplitz matrix. Therefore, the optimized Toeplitz matrices with covariance lags are closer to the sample matrices , since the LR is larger. But are they closer to the true covariance Toeplitz matrix as well? To assess the achieved proximity of the optimized Toeplitz matrix to this true covariance matrix rather than to the sample matrix , let us introduce the "proximity ratio":

Figure 8: fmincon PR Progression with Respect to the True Covariance Matrix ��

The proximity ratio, calculated over 955 iterations for the LR maximization interval in Carathéodory parameters, is followed by 40 iterations of LR maximization over covariance lags. While the LR maximization is performed on the sample matrix, the resulting value is the proximity ratio. Also surprising is the proximity ratio, calculated for the interval with Carathéodory parameters optimization, which remains practically equal to the LR values over this interval.

More surprising are the results of the "proximity ratio" analysis over the second step of the LR maximization, performed over the covariance lags of the optimized Toeplitz matrix. Indeed, the results of this second-step LR maximization led to a degradation of the "proximity ratio", with the losses in PR being practically equal to the LR gain achieved in the second step of LR maximization of covariance lag parameters.

It seems rather natural that all the further gains in LR values, provided by fmincon using the covariance lags for the Toeplitz matrix representation, were lost when the received solution was compared with the original true Toeplitz covariance matrix and not with the sample matrix , as per the likelihood ratio. This property evokes certain doubts in the application of the "second step" LR maximization using the covariance lag parameters for the optimized Toeplitz matrix parameterization. This concern is illustrated by another example of a successful trial (Figure 9), which converged to the solution with LR = 1.8529 ^. 10−1 after 481 iterations.

Figure 9: fmincon LR Progression During the Initial Stage of Optimization in “DoA” Parameters, Followed by Optimization in “Natural” Parameters

The true Toeplitz covariance matrix generated the likelihood ratio LR = 1.494 ^. 10−1. When the true covariance matrix was used to initialize the next fmincon optimization, the LR remained the same, LR = 1.8529 ^. 10−1. Therefore, the 180th initialization was successful, following 179 previous unsuccessful ones, that ended with negative minimal eigenvalues but reached the same Maximum Likelihood value (LR = 1.8529 ^. 10−1) as the fmincon routine when launched from the true covariance matrix, which had LR = 1.494 ^. 10−1. This experiment was then repeated with different random- like initializations of the fmincon routine.

The Carathéodory parameterization did not produce solutions with negative eigenvalues, but in the overwhelming majority of conducted trials, it did not achieve LR values comparable to those from a successful (i.e., with no negative eigenvalues) LR maximization that exploited the covariance lags parameterization. Therefore, for the class of Toeplitz covariance matrices with an unspecified number of equal minimal eigenvalues, the proposed optimization of the "oversampled" (T ≥ N) LR by initial optimization in Carathéodory parameters with m = N − 1 point sources, followed by optimization over the Toeplitz matrix's covariance lags, provides the best convergence rate to the LR maximum, measured by the number of required initializations.

Numerical Likelihood Ratio Maximization for the Models with Limited Signal Subspace Dimension

In the important class of direction of arrival estimation problems, the covariance matrix has a finite dimension of its signal subspace. The number of active sources m (signal subspace dimension) is usually estimated before DOA estimation and is treated as a known parameter. Therefore, for the class of uniform linear arrays (ULAs), the DOA estimation of sources is practically identical to the problem of Hermitian Toeplitz p.d. matrix reconstruction, with the Toeplitz matrix parameterized as in (27) with .

For this type of problem of DOA estimation of sources, the MUSIC solutions are often produced for a sample volume T that is smaller than the matrix dimension N and exceeds the number of sources m:

                                                           m ≤ T < N. (33)

Therefore, the Maximum Likelihood solution should be developed for the "undersampled" training data (T < N) as well. The requirement to expand the maximum likelihood Toeplitz matrix estimation over the class of Toeplitz matrices with a finite signal subspace dimension is an important distinction with respect to the problems addressed in the previous chapter.

Another important distinction with the previous problem is the known limited signal subspace dimension m. Since this type of Toeplitz covariance matrix is very similar to (27):

it is quite natural to use this model instead of (27) with m = N − 1 for LR maximization. Yet, the model (27) with its m = N − 1 describes the full class of p.d. Toeplitz Hermitian matrices, as demonstrated above. Therefore, the reduction of the "full" model (27) to the model (34) with the smaller may affect our ability to calculate the maximum likelihood Toeplitz matrix estimate. The Toeplitz matrix parameterization based on its elements is hardly applicable in this case.

Indeed, the constraints on equality of (N− m ) minimal eigenvalues of the matrix may be as difficult to retain, if not more difficult, than the condition on the positive definiteness of the optimized Toeplitz matrix. Our attempts to use the model (34) with for the MATLAB fmincon optimization failed, since the iterations converged to an LR value that remained smaller than the LR value of the true covariance matrix . A mathematically rigorous explanation of this phenomenon is still missing, but a few negative examples suffice to demonstrate this property.

Note that for the "oversampling" training conditions (T ≥ N), the eigenvalues of the sample covariance matrix are all different with probability one. Therefore, a different approach can be applied to the problem of maximum likelihood Toeplitz matrix estimation with a finite dimension m < N of the matrix signal subspace. For T> N, the sample matrix has all different eigenvalues, and it is possible to look for the maximum likelihood Toeplitz matrix with all different eigenvalues as well.

It is clear that the global maximum likelihood of this solution should exceed the global maximum likelihood for the finite signal subspace dimension. Since such an ML Toeplitz matrix estimate may always be presented as the covariance matrix of (N − 1) sources plus white noise, one can present the received solution in this format and then select the m most powerful sources. The likelihood of this solution must exceed that of the true covariance matrix, which can be checked during Monte Carlo simulations. This approximation can be found by alternating projections over the "full" rank Toeplitz solutions presented in format (25), by the straightforward exclusion of the (N −1− m ) weakest sources.

For applications with the "undersampled" training condition (T < N), the optimization of the Toeplitz matrix should be different. One approach is to find the "diagonally loaded" Toeplitz matrix:

where T(m ) is the Toeplitz matrix of rank m. For this matrix, we may first try to find the Toeplitz p.d. matrix that maximizes the LR ratio modified in [1,16]:

and then find the best approximation of this Toeplitz matrix by a p.d. Toeplitz matrix with a smaller signal subspace dimension m < N. A special investigation is required to determine the optimal m for signal subspace dimension estimation. It is very likely that for applications with a ULA receiver array, better performance may be achieved if the signal subspace dimension m is identified after the maximum likelihood Toeplitz matrix is estimated. While the classical MUSIC may be applied to this ML Toeplitz matrix estimate, in this study, we tested alternating projections to find the p.d. Toeplitz matrix of dimension m for the signal subspace, avoiding the "full" Toeplitz matrix transformation to the "Carathéodory" format (25).

Let us first test the "oversampled" case (T= 85) of DOA estimation by comparing the results of the classical MUSIC applied to the sample matrix to the ML Toeplitz matrix full rank solution.

The ML Toeplitz matrix estimate was obtained by first optimizing the Carathéodory representation (25) with the ultimate number of (N−1) sources, followed by ML optimization of the Toeplitz matrix's covariance lags.

We consider the scenario with m = 3 independent sources acting upon an N=17-element uniform linear array with d/= 0.5. The three sources have DOAs:

The results of the two-step LR optimization are illustrated in Figure 10, which demonstrates the worst estimated DOA error over the three sources. Recall that in the first step, we used the Carathéodory representation, and in the second step, the Vandermonde representation. The analysis of the "proximity ratio" with respect to the true covariance matrix (instead of the sample matrix) is presented on Figure 11. It demonstrates a much higher value of the "proximity" to the true covariance matrix than to the sample matrix in the maximized likelihood ratio.

Figure 10: Maximal “DoAs” Estimation Error PDF

Figure 11: fmincon Proximity Ratio Progression for LR-optimized “DoAs” Reconstructed Toeplitz Covariance Matrix

Even though the final step of LR maximization over the covariance lags of the optimized Toeplitz matrix increased the likelihood ratio from LR = 1.630 ^. 10−1 to LR = 1.863 ^. 10−1, the "proximity ratio" in fact degraded as a result of this optimization from LR = 0. 916, achieved by the Carathéodory parameters optimization, to LR = 0.802 produced by it due to the second-stage LR maximization in covariance lags.

This phenomenon was observed across all conducted trials, suggesting that further LR enhancement beyond the level achieved by the Carathéodory parameterization is not helpful for DOA estimation. Note that for the selected parameters, the classical MUSIC applied to the sample matrix does not resolve the second and third sources in any of the 1,000 trials, while in the proposed Maximum Likelihood Toeplitz matrix estimation, all three sources were resolved in every trial.

In Figure 12, we provide the conventional MUSIC pseudospectrum demonstrating this phenomenon.

Figure 12: 1,000 MUSIC Pseudo-Spectra of the Second and Third (unresolved) Sources of the Original Sample Matrix

Let us continue by presenting the results of maximum likelihood Toeplitz matrix estimation by analyzing the maximum (over the three acting sources) estimation error. Since LR optimization over the covariance lags increased the likelihood ratio but decreased the "proximity" of the solution to the true Toeplitz matrix , let us separately analyze the results of the first stage of LR maximization using the Carathéodory parameterization and the results of the following LR optimization of the covariance lags. As expected, despite the second stage of LR maximization in covariance lags increasing the LR values, the DOA estimation accuracy did not improve.

This result is expected because the "proximity factor" above degraded during the second stage of LR maximization by fmincon, which used the covariance lags as the optimized parameters. The average over all three acting sources RMSE is equal to RMSE = 0.1617o with the standard deviation = 0.1006 for the first-stage LR optimization in Carathéodory parameters, and RMSE = 0.1625o and = 0.1006 for the second stage of LR improvement via covariance lag optimization. One can see that the improvement in the LR value at the second stage of optimization did not translate into improved DOA estimation accuracy.

Obviously, Monte-Carlo simulation results should be compared with the associated Cramér-Rao lower bounds. The Fisher information matrix (FIM) is calculated using the formula [26]:

Recall that the traditional MUSIC, applied to the sample covariance matrix, failed to resolve these close-in sources. In Figure 14 and Figure 15, we introduce 1,000 MUSIC pseudospectra, calculated for Toeplitz matrices with maximized LR, first using the Carathéodory Toeplitz matrix parameterization and then maximizing LR using the covariance lags directly.

Figure 14: 1,000 MUSIC Pseudo-Spectra of the First (top figure) and Second and Third Sources (bottom figure) after LR Maximization in “DoAs”-Only Parameters

Figure 15: 1,000 MUSIC Pseudo-Spectra of the First (top figure) and Second and Third Sources (bottom figure) after Two-Stage LR Maximization in the “DoAs” and “Traditional” Parameters

Despite the second step of LR maximization using covariance lags slightly increasing the optimized LR values, the DOA estimation accuracy in both steps remains largely unchanged. Note also that the restriction of the signal subspace dimension by alternating projections sufficiently increased the dynamic range of the MUSIC pseudo-spectrum, but did not improve the spectral peak positions and associated DOAs estimation errors.

In Figure 16, we compare the MUSIC pseudo-spectra of the original ML Toeplitz matrix with the results of alternating projections that left only three signal subspace eigenvalues above the noise eigenvalue floor in the modified p.d. Toeplitz matrix.

Figure 16: MUSIC Spectrum: Original ML Toeplitz Matrix vs. Alternating Projections

LR Maximization in the ''Undersampled'' Training Conditions

Conclusions and Recommendations

The analysis of the numerical maximization of the likelihood ratio for a p.d. Toeplitz covariance matrix revealed several properties that we believe are important from both theoretical and practical perspectives.

In particular, we demonstrated that numerical LR maximization, utilizing Vandermonde parameterization, applied to random initial Toeplitz matrices, can often yield an inappropriate solution with negative eigenvalues. Yet, in all cases where the MATLAB fmincon routine converges to a positive definite Toeplitz matrix, the process converges to the same solution as if initiated by a true Toeplitz covariance matrix, irrespective of the initial matrix. Our numerous attempts to initiate the fmincon routine with matrices far from the true covariance matrix with LR values < 10−20 (!), nevertheless, finally converged to the same solution as the one initiated by the true covariance matrix.

The obtained calculation results cannot be treated as proof of global optimality. Yet, the substantial volume of optimization results—obtained from a variety of random initial approximations that converged to the same solution—allows for the cautious conjecture that the optimization problem, as formulated here, is convex and that the identified optimum is global.

Naturally, this assertion requires rigorous mathematical proof. Indeed, the fact that—starting from the true covariance matrix— we converge to the very same solution serves as a significant argument in support of such a proof. A further argument in favor of this proof is that this solution yields a likelihood function value that exceeds that of the true covariance matrix and coincides with the value attained when the process is initialized with the true Toeplitz covariance matrix.

Note that the finite calculation accuracy leads to some insignificant variations of the optimum solution, and this effect needs to be closely monitored for larger ULA dimensions (N > 17) and eigenvalue spread (> 10¹⁰).

The proposed two-step optimization starts from the Carathéodory parameterization and, upon convergence, switches to the Vandermonde parameterization. While the maximized LR increases during this second stage of optimization, the proximity ratio relative to the true matrix decreases.

The numerical analysis discussed above was conducted only for N = 17; therefore, the revealed properties should be closely monitored for uniform linear arrays with larger apertures and larger dynamic range, ₁/_N. This particular sequence of Toeplitz matrix parameterizations was shown to be appropriate for the LR maximization of Toeplitz matrices when the Toeplitz matrix has a known a priori number of noise subspace eigenvalues. We demonstrated that the Toeplitz p.d. matrices with a known number of noise subspace eigenvalues cannot be optimized by the fmincon routine.

Specifically, such an LR maximization stops before the LR of the true covariance matrix is reached. For this reason, the problem of LR maximization of the p.d. Toeplitz matrix with a priori known finite signal subspace dimension is proposed to be solved in two steps. First, the problem of maximum likelihood Toeplitz matrix estimation is resolved with no control over the reduction of the signal subspace dimension. In the second stage, for the resulting p.d. Toeplitz matrix with a (globally) maximum-likelihood ratio and distinct eigenvalues, the second problem of noise subspace eigenvalue equalization is addressed using the method of alternating projections. For example, we demonstrated that noise eigenvalue equalization somewhat reduces the maximum likelihood ratio of the final optimized Toeplitz matrices. However, the remaining likelihood ratio still exceeds that of the true Toeplitz covariance matrix.

Moreover, the proposed two-stage version of the likelihood ratio maximization procedure enabled us to significantly expand the method's scope of applicability— encompassing sources situated closer to one another, as well as weaker sources—and to successfully perform their separation and the optimal estimation of their directions of arrival, whereas the standard MUSIC algorithm, when applied to the Hermitian (sample) maximum likelihood matrix, fails to handle this task [28-52].

References

1. Abramovich, Y., Abramovich, V., & Pongsiri, T. (2025). Maximum Likelihood Toeplitz Covariance Matrix Estimation: Is It a Convex Optimization Problem?. Authorea Preprints.
2. Abramovich, Y., Abramovich, V., & Pongsiri, T. (2026). Numerical Techniques for the Maximum Likelihood Toeplitz Covariance Matrix Estimation: Part II. Hermitian Toeplitz Matrices. IEEE Transactions on Aerospace and Electronic Systems, 62, 3971-3990.
3. Burg, J. P., Luenberger, D. G., & Wenger, D. L. (1982). Estimation of structured covariance matrices. Proceedings of the IEEE, 70(9), 963-974.
4. Fuhrmann, D. R. (1991). Application of Toeplitz covariance estimation to adaptive beamforming and detection. IEEE Transactions on Signal Processing, 39(10), 2194-2198.
5. Fuhrmann, D. R., & Barton, T. A. (1990, October). Estimation of block-Toeplitz covariance matrices. In 1990 Conference Record Twenty-Fourth Asilomar Conference on Signals, Systems and Computers, 1990. (Vol. 2, p. 779). IEEE.
6. Robey, F. C. (1990). A covariance matrix modeling approach to adaptive beamforming and detection (Doctoral dissertation, PhD dissertation, Washington Univ., St. Louis, MO).
7. Fuhrmann, D. R., Turmon, M. J., & Miller, M. I. (1988, March). Efficient implementation of the EM algorithm for Toeplitz covariance estimation. In Proceedings of the Twenty-Second Annual Conference on Information Sciences and Systems.
8. Miller, M. I., Fuhrmann, D. R., O’Sullivan, J. A., & Snyder,D. L. (1991). Maximum-likelihood methods for Toeplitz covariance estimation and radar imaging. Advances in Spectrum Analysis and Array Processing, 2, 145-172.
9. Turmon, M. J., Miller, M. I., Snyder, D. L., & O'Sullivan,J. A. (1988, August). Performance evaluation of maximumlikelihood Toeplitz covariance estimates generated using the EM algorithm. In Proc. Fourth ASSP Workshop on Spect. Est. and Modeling (pp. 182-185).
10. Turmon, M. J., & Miller, M. I. (2002). Maximum-likelihood estimation of complex sinusoids and Toeplitz covariances. IEEE transactions on signal processing, 42(5), 1074-1086.
11. Gray, D. A., Anderson, B. D. O., & Sim, P. K. (1987).Estimation of structured covariances with application to array beamforming. Circuits, Systems and Signal Processing, 6(4), 421-447.
12. Vinogradova, J., Couillet, R., & Hachem, W. (2015). Estimation of Toeplitz covariance matrices in large dimensional regime with application to source detection. IEEE Transactions on Signal Processing, 63(18), 4903-4913.
13. Cai, T. T., Ren, Z., & Zhou, H. H. (2013). Optimal rates of convergence for estimating Toeplitz covariance matrices. Probability Theory and Related Fields, 156(1), 101-143.
14. Vallet, P., & Loubaton, P. (2014, May). Toeplitz rectification and DOA estimation with music. In 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) (pp. 2237-2241). IEEE.
15. Abramovich, Y. I., & Spencer, N. K. (2005, May). Expected-likelihood covariance matrix estimation for adaptive detection. In IEEE International Radar Conference, 2005. (pp. 623-628). IEEE.
16. Abramovich, Y. I., & Johnson, B. A. (2007, December). Expected likelihood estimation: Asymptotic properties for" stochastic" complex Gaussian models. In 2007 2nd IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (pp. 33-36). IEEE.
17. Abramovich, Y. I., Arov, D. Z., Kachur, V. G. (1988). Adaptive cancellation filters for stationary interference with a Toeplitz correlation matrix. J. Comm. Tech. Elect., 33(4), 54-61.
18. Aubry, A., Babu, P., De Maio, A., & Rosamilia, M. (2024). Advanced methods for MLE of Toeplitz structured covariance matrices with applications to RADAR problems. IEEE Transactions on Information Theory, 70(12), 9277-9292.
19. Yang, Z., Chen, X., & Wu, X. (2023). A robust and statistically efficient maximum-likelihood method for DOA estimation using sparse linear arrays. IEEE Transactions on Aerospace and Electronic Systems, 59(5), 6798-6812.
20. Cederberg, D. (2024). Toeplitz covariance estimation with applications to MUSIC. Signal Processing, 221, 109506.
21. Muirhead, R. J. (2009). Aspects of multivariate statistical theory. John Wiley & Sons.
22. Abramovich, Y. I., Spencer, N. K., & Gorokhov, A. Y. (2004). Bounds on maximum likelihood ratio-Part I: Application to antenna array detection-estimation with perfect wavefront coherence. IEEE transactions on signal processing, 52(6), 1524-1536.
23. MathWorks. (2024). Problem when using fmincon: The solver starts from a feasible point but converges to an infeasible point.
24. Kolosov, A. B. (1987). Over-the-Horizon Radar. MA, Artech House.
25. Valmarde, J. (2015). Generating symmetric positive definite Toeplitz matrices.
26. Skolnik, M., Nemhauser, G., & Sherman, J. (1964). Dynamic programming applied to unequally spaced arrays. IEEE Transactions on antennas and propagation, 12(1), 35-43.
27. Nagarsenker, B. N., & Pillai, K. C. S. (1973). The distribution of the sphericity test criterion. Journal of Multivariate Analysis, 3(2), 226-235.
28. Horn, R. A., Johnson, C. R. (1985). C. Johnson, Matrix Analysis. Cambridge.
29. Abramovich, Y., Danilov, B.G. (1976). The use of dynamic programming for the synthesis of array antennas with discrete phase shifters. Radio Eng. Elect. Phys., 21(1), 69-74.
30. Paulraj, A., & Kailath, T. (1985, April). Direction of arrival estimation by eigenstructure methods with unknown sensor gain and phase. In ICASSP'85. IEEE International Conference on Acoustics, Speech, and Signal Processing (Vol. 10, pp. 640-643). IEEE.
31. Weiss, A., & Yeredor, A. (2020). Asymptotically optimal blind calibration of uniform linear sensor arrays for narrowband Gaussian signals. IEEE Transactions on Signal Processing, 68, 5322-5333.
32. Weiss, A., & Yeredor, A. (2021, June). Enhanced blind calibration of uniform linear arrays with one-bit quantization by kullback-leibler divergence covariance fitting. In ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) (pp. 4625-4629). IEEE.
33. Chu, M. T. (1993). The stability group of symmetric Toeplitz matrices. Linear algebra and its applications, 185, 119-123.
34. MathWorks. fmincon - Nonlinear Optimization.
35. Chu, M. T., Golub, G. H. (2002). Inverse eigenvalue problems. Oxford Science Publications, 1-71.
36. Shan, T. J., & Kailath, T. (1985). Adaptive beamforming for coherent signals and interference. IEEE Transactions on Acoustics, Speech, and Signal Processing, 33(3), 527-536.
37. Paulraj, A., Reddy, V. U., Shan, T. J., & Kailath, T. (1986, October). Performance analysis of the MUSIC algorithm with spatial smoothing in the presence of coherent sources. In MILCOM 1986-IEEE Military Communications Conference: Communications-Computers: Teamed for the 90's (Vol. 3, pp. 41-45). IEEE.
38. Linebarger, D. A., & Johnson, D. H. (2002). The effect of spatial averaging on spatial correlation matrices in the presence of coherent signals. IEEE transactions on acoustics, speech, and signal processing, 38(5), 880-884. 
39. Indukumar, K. C., & Reddy, V. U. (1992). Anote on redundancy averaging. IEEE Transactions on signal processing, 40(2), 466-469.
40. Shan, T. J., Wax, M., & Kailath, T. (1985). On spatial smoothing for direction-of-arrival estimation of coherent signals. IEEE Transactions on Acoustics, Speech, and Signal Processing, 33(4), 806-811.
41. Doron, M. A., & Weiss, A. J. (2002). Performance analysis of direction finding using lag redundancy averaging. IEEE transactions on signal processing, 41(3), 1386-1391.
42. Turmon, M. J. (1994, June). Cramér-Rao Bounds for Toeplitz Covariance Estimation.
43. Consul, P. C. (1969). The exact distribution of likelihoodcriterion of different hypotheses. In Multivariate Analysis II,P.R. Krishnain, Ed., Academic, 171-181.
44. Gradshteyn, L. S., Ryzhik, L. M. (2000). Tables of integrals, series, and products. Sixth Edition, Academic Press.
45. Mestre, X. (2008). Improved estimation of eigenvalues and eigenvectors of covariance matrices using their sample estimates. IEEE Transactions on Information Theory, 54(11), 5113-5129.
46. MathWorks. linprog – solve linear programming problems.
47. Wu, W. B., & Pourahmadi, M. (2009). Banding sample autocovariance matrices of stationary processes. Statistica Sinica, 1755-1768.
48. Chun, J., & Kailath, T. (1989). A constructive proof of the Gohberg-Semencul formula. Linear Algebra and its Applications, 121, 475-489.
49. Evans, J. E. (1981). High resolution angular spectrum estimation technique for terrain scattering analysis and angle of arrival estimation. In 1st IEEE ASSP Workshop Spectral Estimat., McMaster Univ., Hamilton, Ont., Canada, 1981 (pp. 134-139).
50. Greenbaum, A., Li, R. C., & Overton, M. L. (2020). First-order perturbation theory for eigenvalues and eigenvectors. SIAM review, 62(2), 463-482.
51. Abramovich, Y., Abramovich, V., Pongsiri, T. (2025). Numerical Techniques for the Maximum Likelihood Toeplitz Covariance Matrix Estimation: Part I. Symmetric Toeplitz Matrices. J Electr Comput Innov, 2(2), 01-21.
52. Abramovich, Y. I., & Johnson, B. A. (2008). GLRT-based detection-estimation for undersampled training conditions. IEEE Transactions on Signal Processing, 56(8), 3600-3612.