# An Efficient Constrained Learning Algorithm for Stable 2D IIR Filter Factorization

Authors: Nicholas Ampazis; Stavros J. Perantonis
Journal: Advances in Artificial Neural Systems
 (2013)

---

## Abstract

A constrained neural network optimization algorithm is presented for factorizing simultaneously the numerator and denominator polynomials of the transfer functions of 2-D IIR filters. The method minimizes a cost function based on the frequency response of the filters, along with simultaneous satisfaction of appropriate constraints, so that factorization is facilitated and the stability of the resulting filter is respected.


---

## Body

## 1. Introduction

Factorization of 2D polynomials is an important problem in the design of IIR filters in two-dimensional signal processing because a factorized filter transfer function can be efficiently implemented in cascade form. However, the fundamental theorem of Algebra, concerning polynomial factorization, cannot be extended to the two-dimensional case; that is, a bivariate polynomialF(z1,z2) cannot, in general, be factored as a product of lower-order polynomials. The 2D IIR cascade structure may be considered as an attractive circuit design alternative, because of its relative insensitivity to coefficient quantization, the requirement of smaller number of arithmetic operations for a given filter size, and finally, because issues related to the stability of the filters are easier to deal with for filters with a smaller number of denominator coefficients [1]. Hence, an efficient method for the factorization of 2D IIR filters would be most beneficial.Many methods have been proposed for dealing with the two-dimensional polynomial factorization problem [2–6]. Most of those methods, however, adopt a conventional numerical method of roots finding (e.g., Laguerre’s, Newton-Raphson’s, Jenkins-Traub’s methods, etc.), in which successive approximations to the roots are obtained. In addition, almost all numerical methods can only find the roots one by one, that is, by deflation method, where the next root is obtained by the deflated polynomial after the former root is found. This means that numerical root-finding methods are inherently sequential.In an earlier work we have proposed a neural network-based Constrained Learning Algorithm (CLA) for factorizing 2D polynomials using constrained optimization techniques [7]. Using that approach we were able to obtain exact solutions for factorable polynomials and excellent approximate solutions for nonfactorable polynomials. The technique offers the advantage of incorporating a priori information about the relations between the coefficients of the original polynomial and the coefficients of the desired factor polynomials. By incorporating additional stability constraints into the formalism, the stability (in the sense of Huang [8]) of a transfer function with the resulting factorized polynomial in the denominator (of up to second order) is also made possible. Our CLA method has been explored further in [9] by utilizing root moments, in [10] with an adaptive neural network architecture, and in [11] for handling complex root cases.Despite its success in 2D polynomial factorization, CLA suffers from disadvantages which can limit its effectiveness for its application to IIR filter factorization. CLA’s cost function is not directly related to a proximity measure between the frequency response characteristics of the filters, but is based on the Mean Square Error (MSE) criterion between the output of the original polynomial and the output of the resulting factorable polynomial for various real values of|z1|<1 and |z2|<1. Thus, CLA requires separate factorizations for the numerator and denominator polynomials of the IIR filter, while better solutions, in terms of the proximity of filter characteristics, could be achieved by simultaneous variation of all available parameters (numerator and denominator coefficients). In this paper we apply constrained optimization techniques to the more general problem of implementing a stable cascade IIR filter, given the specifications of an original 2D filter. This is achieved by suitably modifying our previous formalism as it is explained in the following section.


## 2. Derivation of the Algorithm

Consider the polynomial(1)f(z1,z2)=∑i=0NA‍∑j=0NA‍aijz1iz2jwithNA even and a00=1. For the above polynomial, we seek to achieve an exact or approximate factorization of the form
(2)f(z1,z2)≈∏i=1,2f(i)(z1,z2),
where
(3)f(i)(z1,z2)=∑j=0MA‍∑k=0MA‍wjk(i)z1jz2k
with MA=NA/2.To this end, we can utilize a feedforward neural network architecture in order to learn the appropriate weightswjk(i) as we proposed in [7]. For example, Figure 1 shows the architecture of such a network for the case where NA=2. The units in the hidden layer have a logarithmic activation function ln, and the output unit is linear. This architecture has been preferred to the ∑‍-∏‍ architecture employed in [12], because it produces smoother cost function landscapes avoiding deep valleys and thus facilitates learning. As we showed in [7], this neural network can be trained using a Constrained Learning Algorithm (CLA) that achieves minimization of the MSE criterion along with simultaneous satisfaction of multiple equality and inequality constraints between the polynomial coefficients. Using that method, we are able to obtain very good approximate solutions for nonfactorable polynomials.Figure 1
Neural network architecture used to solve the factorization problem.Now consider the 2D IIR filter transfer function(4)H(z1,z2)=B(z1,z2)A(z1,z2),
where
(5)B(z1,z2)=∑i=0NB‍∑j=0NB‍bijz1iz2j,A(z1,z2)=∑i=0NA‍∑j=0NA‍aijz1iz2j
with NA, NB even and b00=a00=1.For the above polynomials, we seek to achieve an exact or approximate factorization of the form(6)B(z1,z2)≈∏i=1,2B(i)(z1,z2),
where(7)B(i)(z1,z2)=∑j=0MB‍∑k=0MB‍ujk(i)z1jz2k,(8)A(z1,z2)≈∏i=1,2A(i)(z1,z2),
where
(9)A(i)(z1,z2)=∑j=0MA‍∑k=0MA‍vjk(i)z1jz2k
with MA=NA/2, and MB=NB/2 so that the transfer function of (4) can be approximated by
(10)H^(z1,z2)=∏i=1,2‍B(i)(z1,z2)∏i=1,2‍A(i)(z1,z2).Unfortunately, if it were to apply CLA in this problem we would have to derive two separate factorizations for the numerator and denominator polynomials. Such an optimization scheme, however, is not directly related to a proximity measure between the original and the factorized filter frequency response characteristics. In terms of the proximity of filter characteristics, obviously a better approximation can be achieved by a simultaneous variation of both the numerator and denominator coefficients.In this case, we should ask: in what sense areH(z1,z2) and H^(z1,z2) similar? The answer to this question is that we would at least expect the amplitude spectra of H(z1,z2) and H^(z1,z2) to be roughly equal so that the specifications of the original and of the factorized filter match as close as possible. Hence, the main purpose of the factorization algorithm is to minimize an MSE cost function based on the magnitude of the frequency responses of the original and the resulting filter:
(11)E=12∑θ1=-ππ‍∑θ2=-ππ‍(T-O)2,T=ln(|H(ejθ1,ejθ2)|),O=ln(|H^(ejθ1,ejθ2)|).
The double summation is carried over a regular grid of points. Logarithmic functions are used also in this case in order to suppress deep valleys and thus facilitate convergence. However, E is still a highly nonlinear function of the polynomial coefficients, and therefore gradient-based techniques can run into problems in the presence of valleys or ridges in the cost function landscape. For such cases it has proved beneficial to incorporate into the adaptation rule additional a priori problem-specific information. This often facilitates the convergence of the algorithm and leads to better solutions.For the problem at hand, additional information is available in the form of constraints between the coefficients of the given numerator and denominator polynomials and the coefficients of the corresponding desired factor polynomials. More explicitly, if we assume thatB(z1,z2) and A(z1,z2) are factorable, then these constraints can be expressed as follows:
(12)Φj+(NB+1)iu=bij-∑l=1i‍∑m=1j‍ulm(1)ui-l,j-m(2)=0,0≤i≤NB,0≤j≤NB,Φj+(NA+1)iv=aij-∑l=1i‍∑m=1j‍vlm(1)vi-l,j-m(2)=0,0≤i≤NA,0≤j≤NA.
Thus, for the general case in which B(z1,z2) and A(z1,z2) are nonfactorable, the objective of the adaptation process is to reach a minimum of the cost function of (11) with respect to the variables ujk(i) and vjk(i), which satisfies as best as possible the constraint Φ=0, whereΦT=(Φu,Φv) with
(13)Φu=(Φj+(NB+1)iu,0≤i≤NB,0≤j≤NB),Φv=(Φj+(NA+1)iv,0≤i≤NA,0≤j≤NA).The strategy which will be followed in order to solve this problem is similar to the method we proposed in [7], but with appropriate modifications in order to take care of the above constraints, as well as the new cost function that we introduced in (11).At each iteration of our algorithm, the vector of unknown coefficientsw=(u,v)T is changed by dw, so that
(14)dwTdw=(δP)2,
where δP is a small constant. Thus, at each iteration, the search for an optimum new point in thewspace is restricted to a small hypersphere centered at the point defined by the current vector w. If δP is small enough, the changes to E and toΦinduced by changes in the coefficients can be approximated by the first differentials dE and dΦ. Also, at each iteration, we seek to achieve the maximum possible change in |dE|, so that (14) is respected and the change dΦin Φis equal to a predetermined vector quantity δQ, designed to bringΦcloser to its target (zero).Maximization ofdE can be carried out analytically by introducing suitable Lagrange multipliers. Thus, a vectorνof Lagrange multipliers is introduced for the constraints in (12) and a multiplier μ is introduced for (14).Consider a functionε, whose differential is defined as follows:
(15)dε=dE+(δQT-dΦT)ν+μ[(δP)2-dwTdw].
On evaluating the differentials involved in the right-hand side, we readily obtain
(16)dε=dwTJ+(δQT-dwTFT)ν+μ[(δP)2-dwTdw],where(i)
J is a vector with elements Ji=∂E/∂wi;(ii)
Fis a matrix given by F=[Fu00Fv];(iii)
Fu andFvare matrices with elements Fmiu=∂Φmu/∂ui and Fmiv=∂Φmv/∂vi, respectively. To minimize dε (maximize its magnitude) at each iteration, we demand that
(17)d2ε=J-FTν-2μdw=0,d3ε=-2μ>0.From (17) it immediately follows that
(18)dw=J2μ-FTν2μ.
Equation (18) is the update rule for the coefficients of the factor polynomials, but μ andνmust still be evaluated in terms of known quantities.Noting thatdΦ=δQ=Fdw,we can multiply both sides of (18) byF and solve for νto obtain
(19)ν=-2μIFF-1δQ+IFF-1IJF,whereIFF=FFT,IJF=FJ.
To evaluate μ, we substitute (18) and (19) into (14) and obtain
(20)μ=-12[IJJ-IJFTIFF-1IJF(δP)2-δQTIFF-1δQ]1/2,whereIJJ=JTJ.
Note that due to the block-diagonal form of matrixF,it is possible to avoid the tedious evaluation of IFF-1 by rewriting the quantities IJJ, IFF, and IJF as follows:  (i)
IJJ=(Ju)TJu+(Jv)TJv, where Jiu=∂E/∂ui and Jiv=∂E/∂vi,(ii)
IFF=IFFu+IFFv=Fu(Fu)T+Fv(Fv)T,(iii)
IJF=IJFu+IJFv=FuJu+FvJv.Hence, (19) and (20) can be rewritten as
(21)ν=-2μ(IFFu)-1δQu-2μ(IFFv)-1δQv+(IFFu)-1IJFu+(IFFv)-1IJFv,μ=-12[IJJ-(IJFu)T(IFFu)-1IJFu-(IJFv)T(IFFv)-1IJFv(δP)2-(δQu)T(IFFu)-1δQu-(δQv)T(IFFv)-1δQv]1/2,
since the matrices (IFFu)-1and (IFFv)-1can be computed more easily than the composite matrixIFF-1.We next discuss the choice ofδQ. It is natural to choose the elements δQm of δQproportional to Φm with a negative constant of proportionality -k, k>0, so that the constraints move towards zero at an exponential rate. Then, the bound k≤δP(ΦTIFF-1Φ)-1/2 set on the value of k by (20) forces us to choose k adaptively. The simplest choice for k will be used; namely, k=[ξ(δP)2/((Φu)T(IFFu)-1Φu + (Φv)T(IFFv)-1Φv)]1/2, where 0<ξ<1. Thus, the final update rule has only two free parameters, namely, δP and ξ.


## 3. Stability Constraints

For the filterH(z1,z2) defined in (4), it is reasonable to assume that for most practical purposes the numerator polynomial B(z1,z2) does not affect the system stability [1]. In the important special case of a second-order denominator polynomial A(z1,z2)(NB=2), the stability conditions according to Huang [8] for the factor polynomials can be written as Xk<1,k=1,…,6, as shown in Table 1. These conditions can be incorporated in the formalism as follows: at a certain iteration of the algorithm, let a number M of the stability constraints be violated, so that Xkr≥1,r=1,…,M. We can augment the vectorΦv of the previous section by adding extra components Φ[j+(NA+1)i]+rv=Xkr,r=1,…,M. By following the formalism of Section 2 using the augmented vectorΦv, the quantities Xkr are driven towards the stability region. Note that the dimensionality of the vectorΦvis adaptive, as only the violated constraints are included in each iteration. Iterative application of this scheme can lead to a final stable solution.Table 1
Stability constraints for the denominator polynomial.
X
1
=
(
v
10
(
1
)
)
2
X
2
=
(
v
10
(
2
)
)
2
X
3
=
(
v
01
(
1
)
)
2
+
(
v
11
(
1
)
)
2
+
2
v
01
(
1
)
v
11
(
1
)
-
(
v
10
(
1
)
)
2
-
2
v
10
(
1
)
X
4
=
(
v
01
(
1
)
)
2
+
(
v
11
(
1
)
)
2
-
2
v
01
(
1
)
v
11
(
1
)
-
(
v
10
(
1
)
)
2
+
2
v
10
(
1
)
X
5
=
(
v
01
(
2
)
)
2
+
(
v
11
(
2
)
)
2
+
2
v
01
(
2
)
v
11
(
2
)
-
(
v
10
(
2
)
)
2
-
2
v
10
(
2
)
X
6
=
(
v
01
(
2
)
)
2
+
(
v
11
(
2
)
)
2
-
2
v
01
(
2
)
v
11
(
2
)
-
(
v
10
(
2
)
)
2
+
2
v
10
(
2
)


## 4. Experimental Results

We seek an approximate factorization of the second-order 2D IIR filter whose numerator and denominator polynomial coefficients are shown in Tables2 and 3, respectively. Note that both polynomials are nonfactorable.Table 2
Factorization results for the numerator polynomial.
Original non-factorable
Proposed method
CLA
Numerator polynomial coefficient matrixB
[
1
0.4000
0.0300
0.2500
0.4500
0.0700
0.0125
0.0350
0.0100
]
[
1
0.4001
0.0384
0.2501
0.4649
0.0809
0.0078
0.0451
0.0427
]
[
1
0.4000
0.0276
0.2500
0.4418
0.0602
0.0156
0.0484
0.0318
]
1st factor polynomial coefficient matrix
—
[
1
0.1598
0.0368
0.1678
]
[
1
0.0888
0.1212
0.1141
]
2nd factor polynomial coefficient matrix
—
[
1
0.2403
0.2133
0.2542
]
[
1
0.3111
0.1287
0.2785
]Table 3
Factorization results for the denominator polynomial.
Original non-factorable
Proposed method
CLA
Denominator polynomial  coefficient matrixA
[
1
-
0.9636
0.2321
-
0.9636
0.8360
-
0.1791
0.2321
-
0.1791
0.0445
]
[
1
-
0.9638
0.2315
-
0.9638
0.8910
-
0.2030
0.2320
-
0.2065
0.0438
]
[
1
-
0.9635
0.2320
-
0.9635
0.8332
-
0.1778
0.2317
-
0.1767
0.0335
]
1st factor polynomial coefficient matrix
—
[
1
-
0.5087
-
0.4649
0.2518
]
[
1
-
0.4902
-
0.4616
0.1625
]
2nd factor polynomial coefficient matrix
—
[
1
-
0.4551
-
0.4989
0.1739
]
[
1
-
0.4733
-
0.5019
0.2062
]Figure2(a) shows the amplitude spectrum of the original low-pass IIR filter, and Figure 3 shows the error curve of the unconstrained gradient descent method applied to (11) using a regular grid of 322 points with an adaptation gain of 10-6. From the latest figure we can see that after an initial error reduction, an undesirable flat minimum is reached followed by unbounded oscillations caused by a transition into a region of the parameter space in which the denominator polynomial is unstable. The amplitude spectrum of the filter corresponding to the factorization result achieved at the flat minimum is shown in Figure 2(b) and is clearly unacceptable.Amplitude spectrum of (a) original filter. (b) Filter obtained with gradient descent. (c) Filter obtained with CLA. (d) Filter obtained with the proposed method.
(a)
(b)
(c)
(d)Figure 3
Error plot of gradient descent method.Using the constraintsΦ1u, Φ3u, Φ4u, Φ1v, Φ3v, and Φ4v and all stability constraints of Section 3, we have applied both CLA, as described in [7], and the method proposed in the present paper using a regular grid of 322 training points. Figure 2(c) shows the amplitude spectrum of the filter obtained from the original CLA approach (with δP=0.01 and ξ=0.99) applied sequentially to the numerator and denominator polynomials. Obviously CLA outperforms gradient descent in all aspects. However, there are still small (but not negligible) differences between the amplitude spectra shown in Figures 2(a) and 2(c).Figure2(d) shows the amplitude spectrum of the filter obtained by the proposed method, with δP=0.07 and ξ=0.99. Clearly, there is very good agreement between the amplitude spectra of Figures 2(a) and 2(d). In order to visualize this better, in Figure 4, we present the differences in the amplitude spectra between the original filter and each of the filters obtained by gradient descent (Figure 4(a)), CLA (Figure 4(b)), and the proposed method (Figure 4(c)). From this figure it is clear that, of all the methods, the one proposed in the present paper achieves the maximum proximity to the original amplitude spectrum, leading to the design of a cascade filter that meets as closely as possible the specifications of the original nonfactorable filter.Differences in the amplitude spectra between the original filter and the resulting factorable filters: (a) Gradient descent. (b) CLA. (c) Proposed method.
(a)
(b)
(c)The factorization results are summarized in Tables2 and 3. Note that the separate factorization of the CLA approach produces coefficients that are closer to those of both the numerator and denominator polynomials of the original filter, while the proposed method achieves maximum similarity between the frequency responses in expense of a smaller coefficient similarity, distributed along both the numerator and denominator coefficients.


## 5. Conclusions

In this paper, an efficient algorithm was proposed for factorizing simultaneously the numerator and denominator polynomials of the transfer functions of 2D IIR filters. The method is useful for designing IIR filters in cascade form for more efficient hardware implementation. By introducing a cost function based on the difference of frequency responses between the original and cascade filters and using a flexible framework for incorporating additional information in the form of constraints, we obtained optimal stable solutions for nonfactorable filters.



---

*Source: 292567-2013-02-24.xml*