A Control Theoretical Approach to Spectral Factorization is Unstable

(1)

https://doi.org/10.1007/s00034-021-01817-3 S H O R T P A P E R

A Control Theoretical Approach to Spectral Factorization is Unstable

Michał Meller¹ ·Adam Lasota¹

Received: 13 April 2021 / Revised: 3 August 2021 / Accepted: 5 August 2021 / Published online: 30 August 2021

Abstract

Local stability analysis of a recently proposed recursive feedback-based approach to spectral factorization is performed. The method is found not to give stability guarantees. Interestingly enough, its global behavior often allows one to obtain reasonable approximations of spectral factorizations if a suitable stopping criterion is employed.

Keywords Spectral factorization·Stability analysis·Associated ordinary differential equation method·Recursive algorithms

1 Introduction

Spectral factorization of polynomials is a classical problem with important applications in control [2], optimal filtering [13,14], and design of invertible filters, among others. In the univariate case, one can summarize the problem of spectral factorization succinctly as follows. LetV(z⁻¹)=v0+v1z⁻¹+ · · · +vNz⁻^Ndenote a polynomial (with real coefficients) of the variablez⁻¹. Denote byV(z)the same polynomial in the variable z, rather thanz⁻¹,V(z)=v0+v1z+ · · · +vNz^N. Note thatz0is a root ofV(z⁻¹) if and only if 1/z0is a root ofV(z). A polynomial A(z⁻¹)is said to be a spectral factorization of the Laurent polynomial V(z⁻¹)V(z)if the following equality takes place

A(z⁻¹)A(z)=V(z⁻¹)V(z) . (1)

B

Michał Meller

michal.meller@pitradwar.com ; michal.meller@eti.pg.edu.pl Adam Lasota

adam.lasota@o2.pl

1 Department of Automatic Control, Faculty of Electronics, Telecommunications and Computer Science, Gda´nsk University of Technology, Narutowicza 11/12, 80-233 Gda´nsk, Poland

(2)

Technically, for a polynomialV(z⁻¹)of degreeN there are up to 2^Nspectral factorizations possible, because replacing any root ofA(z⁻¹)with its reciprocal from A(z) also results in a valid spectral factorization.¹This ambiguity is typically removed by requesting a minimum-phase factorization, i.e., a spectral factorization with all roots inside the unit circle on theZplane.

Minimum-phase spectral factorization of small polynomials is easy nowadays because it can be accomplished simply by numerically computing roots ofV(z⁻¹) and “reflecting” those outside the unit circle into its inside, z0 ← 1/z0. For high- degree polynomials, spectral factorization remains a challenging problem because of high computation complexity and numerical sensitivity of solving for roots directly.

Multiple methods for solving this problem, as well as for solving the more general problem of factoring multivariate rational spectral densities, were developed and the reader is referred to, e.g., [1,3,4,6,11,12] for more details and examples.

Unfortunately, majority of the existing approaches are not suitable for embedded applications that require spectral factorization to be performed online in real time, such as the iterative learning active noise control [7]. A noteworthy exception to this rule is a spectral factorization scheme inspired by control theory proposed in [10]. The method is remarkably simple and can be implemented even on the poorest of microcontrollers.

The spectral factorization is obtained in a process of iterative refinement that includes negative feedback.

Even though it is very well known in the control theory that feedback alone does not guarantee stability, no stability analysis of the algorithm was performed. In fact, our experience with the method showed that it often diverges. In this paper, we show using an analytic approach that the method, in general, does not guarantee stability of the convergence irrespective of the value of the step-size parameter used.

2 Stability Analysis 2.1 Preliminaries

In this paper, we will consider only the univariate case, as doing so makes discussion simpler. To quickly demonstrate the instability in the multivariate case [9] on the basis of the univariate case, it is sufficient to consider multivariate polynomials whose coefficients are multiples of the eye matrix.

LetX(z⁻¹)be a polynomial or a finite Laurent series and denote by [X]₋a column vector that contains only those coefficients of X(z⁻¹)that have nonnegative powers ofz⁻¹, sorted in the ascending order (or, equivalently, nonpositive powers ofzsorted in the descending order). For example, ifX(z⁻¹)= −3z²+2z+1+2z⁻¹−3z⁻², then

X(z⁻¹)

−= [1 2 −3]^T. LetL(z⁻¹)=V(z⁻¹)V(z)be a Laurent polynomial whose spectral factorization is sought. Using the introduced “square brace” notation, one can express the method proposed in [10] as follows

Wk+1(z⁻¹)

−=

Wk(z⁻¹)

−+ηk , (2)

1 If one treatsA(z⁻¹)and−A(z⁻¹)as different factorizations, the number grows to 2^N⁺¹.

(3)

where Wk(z⁻¹)is the, iteratively refined, estimate of the spectral factorization at iterationk,

k =

L(z⁻¹)−Wk(z⁻¹)Wk(z)

− (3)

is the residual vector, andη >0 is a small step-size parameter.

2.2 Stability Analysis

Let A(z⁻¹) be a spectral factorization of L(z⁻¹). It is straightforward to see that A(z⁻¹)is a stationary point of (2), because whenWk(z⁻¹)= A(z⁻¹)one obtains that k =0andWk+1(z⁻¹)=Wk(z⁻¹). Overall, the algorithm has up to 2^N⁺¹stationary points.

We will analyze behavior of (2) in a neighborhood of A(z⁻¹), i.e., for Wk(z⁻¹) close toA(z⁻¹). Set

ΔWk(z⁻¹)=Wk(z⁻¹)−A(z⁻¹) . Then, it is straightforward to obtain that

ΔWk+1(z⁻¹)

−=

ΔWk(z⁻¹)

−+ηk (4)

and

k = −

A(z⁻¹)ΔWk(z)+ΔWk(z⁻¹)A(z)+ΔWk(z⁻¹)ΔWk(z)

− . (5) When expressed in this form, a stationary point of (4) is in the origin.

One can analyze behavior of discrete-time algorithms in the form (4) using the associated ordinary differential equation (ODE) method [5,8]. The ODE associated with (4) has the form

dΔWs(z⁻¹)

ds =s , (6)

wheresis a continuous time variable that corresponds to the rescaled original discrete- time,s∼ηk. For small values of the step-size parameterη, trajectories obtained by solving (6) and (4) converge—see [5,8] for more details.

ForΔWs(z⁻¹)sufficiently close to zero, one may neglect second-order terms in (5), which leads to

s −

A(z⁻¹)ΔWs(z)+ΔWs(z⁻¹)A(z)

− . (7)

(4)

Furthermore, it is straightforward to check thats can be expressed in the form s =M

ΔWs(z⁻¹)

− , (8)

where

M = −(M1+M2) (9) andM1,M2are the Toeplitz and the Hankel matrices of the following forms

M1=

⎡

⎢⎢

⎣

a0 a1 . . . aN

0 a0 . . .aN−1

... ... ... ...

0 0 . . . a0

⎤

⎥⎥

⎦

M2=

⎡

⎢⎢

⎣

a0 a1 . . .aN−1aN

a1 a2 . . . aN 0 ... ... ... ... 0 aN 0 0 . . . 0

⎤

⎥⎥

⎦ . (10)

The associated ODE is locally stable if and only if all eigenvalues ofMhave strictly negative real parts. In this case, (2) is locally stable for sufficiently small gainη. On the other hand, ifM has at least one eigenvalue that has nonnegative real part, (2) is locally unstable around A(z⁻¹)no matter how smallηis [5,8].

Note that, since there could exist more than one spectral factorization of a given L(z⁻¹), to understand the algorithm’s global behavior better, one should carry out the above stability analysis for all possible polynomials A(z⁻¹)that satisfy L(z⁻¹) = A(z⁻¹)A(z). If none of these stationary points is found locally stable, then it can be concluded the algorithm cannot possibly converge to a spectral factorization of L(z⁻¹). If, on the other hand, one or more stationary points are found to be locally stable, convergence is possible, but can only be guaranteed locally.

2.3 Computer Simulations

Analyzing eigenvalues of (9), one can see that the algorithm is guaranteed to be locally stable for positive scalars (in which case it performs a recursive computation of a square root). In case of a first-order polynomial A(z⁻¹), it is straightforward to show that the algorithm is locally stable if and only if A(z⁻¹)is a minimum-phase polynomial.

Establishing analytic conditions of stability for polynomials with degree higher than 1 is difficult, if at all possible, but it relatively easy to find counterexamples that show that the algorithm can be unstable even for minimum-phase polynomials. Consider as an example the polynomialV(z⁻¹)=1+1.7826z⁻¹+0.7944z⁻²and setL(z⁻¹)= V(z⁻¹)V(z). For the stationary pointA(z⁻¹)=V(z⁻¹), the matrixMhas eigenvalues

−4.8435, 0.0246±0.0546j, i.e., our analysis predicts that the method will be locally

(5)

0 2 4 6 8 10 x 10⁴ 1

2 3 4 5x 10⁻³

k

Errornorm

Fig. 1 Error trajectory for the example in the text

unstable. We set the initial conditions toW0(z⁻¹)=1.0007+1.7842z⁻¹+0.7937z⁻², i.e., very close toA(z⁻¹)and let the algorithm (2) run for 100,000 iterations using a very small step-sizeη=2.5·10⁻⁵. The trajectory of the error, defined as

ΔWk(z⁻¹)

−², shown in Fig.1is consistent with the analytic predictions. Moreover, all remaining spectral factorizations ofL(z⁻¹)are unstable stationary points, which means that the algorithm simply cannot converge to a factorization ofL(z⁻¹).

For the case of Example 2 from [10], whereA(z⁻¹)=1+0.4573z⁻²+0.4077z⁻⁴+ 0.3101z⁻⁶+0.11831z⁻⁸−0.2587z⁻¹⁰(note that, in [10] there are several typographic errors in the signs of the polynomial’s coefficients), and the algorithm converges, our method confirms that the stationary point is indeed stable.

We validated our analysis extensively using Monte Carlo simulations. In each experiment, a minimum-phase polynomialA(z⁻¹)with degree between 4 and 20 was generated randomly. The algorithm (2) was initialized by adding small random pertur- bations to coefficients ofA(z⁻¹)(Gaussian distributed, zero-mean, standard deviation equal 10⁻³) and ran for 10,000–100,000 iterations. The behavior of the resultant error trajectory was compared with an analytic prediction regarding the algorithm’s stability obtained from (9). In several thousand such experiments, the results of ODE analysis always agreed with simulations.

3 Global Behavior

Oddly enough, the global behavior of the feedback-based method is somewhat more benign. Very often the algorithm initially manages to converge to a close approximation of a spectral factorization. The instability mechanism shown in Sect.2dominates at later stages and the process eventually diverges. Therefore, monitoring the growth of k²and stopping the algorithm when it starts to grow may allow one to obtain an approximation of spectral factorizations with accuracy that may be sufficient for some applications. Nevertheless, there are simply no guarantees in the algorithm that such behavior will take place and that such a heuristic modification will work reliably.

(6)

4 Conclusions

The control theoretic spectral factorization method from [10] was analyzed using the method of associated ordinary differential equations and it was shown that, in general, it does not give stability guarantees.

Data Availibility The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visithttp://creativecommons.org/licenses/by/4.0/.

References

1. B. Anderson, K. Hitz, N. Diem, Recursive algorithm for spectral factorization. IEEE Trans. Circuits Syst.21(6), 742–750 (1974)

2. K.J. Åström, B. Wittenmark,Adaptive Control(Addison-Wesley Longman Publishing Co., Inc, Boston, 1994)

3. G. Baggio, A. Ferrante, On minimal spectral factors with zeroes and poles lying on prescribed regions.

IEEE Trans. Autom. Control61(8), 2251–2255 (2016)

4. G. Baggio, A. Ferrante, On the factorization of rational discrete-time spectral densities. IEEE Trans.

Autom. Control61(4), 969–981 (2016)

5. A. Benveniste, M. Métivier, P. Priouret,Adaptive Algorithms and Stochastic Approximations(Springer- Verlag, 1990)

6. G. Janashia, E. Lagvilava, L. Ephremidze, A new method of matrix spectral factorization. IEEE Trans.

Inf. Theory57(4), 2318–2326 (2011)

7. A. Lasota, M. Meller, Iterative learning approach to active noise control of highly autocorrelated signals with applications to machinery noise. IET Signal Process.14, 560–568 (2020)

8. L. Ljung, Analysis of recursive stochastic algorithms. IEEE Trans. Autom. Control22(4), 551–575 (1977).https://doi.org/10.1109/TAC.1977.1101561

9. T. Moir, Control theoretical approach to multivariable spectral factorisation problem. Electron. Lett.

45(1), 1215–1216 (2009)

10. T. Moir, A control theoretical approach to the polynomial spectral-factorization problem. Circuits Syst.

Signal Process.30, 987–998 (2011).https://doi.org/10.1007/s00034-011-9270-4

11. J. Rissanen, Algorithms for triangular decomposition of block hankel and toeplitz matrices with appli- cation to factoring positive matrix polynomials. Math. Comput.27(121), 147–154 (1973)

12. A.H. Sayed, T. Kailath, A survey of spectral factorization methods. Numer. Linear Algebra Appl.

8(6–7), 467–496 (2001).https://doi.org/10.1002/nla.250

13. H.L. van Trees, K.L. Bell, Z. Tian,Detection, Estimation and Modulation Theory (Detection, Estima- tion, and Filtering Theory, Part I(Wiley, 2013)

14. N. Wiener,Extrapolation, Interpolation and Smoothing of Stationary Time Series(Wiley, New York, 1949)

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.