Jacobian Computation in Least Squares IIR Filter Design

Question

Long time lurker and first time poster - but unfortunately I haven't had any joy untangling this on my own.

I've been studying Mathias Lang's thesis, Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Response, with particular interest in the method for least squares design of stable IIR filters.

However, I'm struggling to understand the Jacobian calculation stage within the mpiir_l2 routine.

Looking through the code, I'm comfortable with the calculations up to the section immediately prior, where the computation of H and the complex response error is completed. After this, however, I fail to follow the Jacobian derivation.

I've included an excerpt as far as the line where my understanding runs out.

function [b,a,l2error] = mpiir_l2(M,N,om,D,W,r,a0)
% MPIIR_L2: [b,a,l2error] = mpiir_l2(M,N,om,D,W,r,a0)
% Least Squares Digital IIR Filter Design with Arbitrary Magnitude
% and Phase Responses and Specified Maximum Pole Radius
%
% OUTPUT:
% b         numerator polynomial coefficients
% a         denominator polynomial coefficients
% l2error   approximation error
%
% INPUT:
% M         order of numerator polynomial
% N         order of denominator polynomial
% om        frequency grid, 0<=om<=pi
% D         complex desired frequency response on the grid om
% W         positive weighting function on the grid om
% r         maximum pole radius
% a0        initial denominator guess; optional
%
% EXAMPLE:
% Lowpass filter with approximately linear passband phase
% (passband group delay = 19 samples, maximum pole radius = 0.97)
% om=pi*[linspace(0,.2,20),linspace(.25,1,75)];
% D=[exp(-j*om(1:20)*19),zeros(1,75)];
% W=[ones(1,20),100*ones(1,75)];
% [b,a,e]=mpiir_l2(16,6,om,D,W,.97);
%
% Author: Mathias C. Lang, Vienna University of Technology, Oct. 98
% uses: lslevin, levin, update, locmax,
%       qp (Optimization Toolbox)

om=om(:); D=D(:); W=W(:); srW = sqrt(W);
EM = exp(-j*om*(0:max([M,N]))); tol = 1e-4; alpha = 0.5;
MAXIT = 100; loopcnt = 0;

fprintf('\n   ERROR  MAX. RADIUS   STEP      SLOPE\n');
fprintf('------------------------------------------\n');

% initial solution
ini=0;
if nargin == 7,
   a=a0(:);
   if length(a)~=N+1, ini=1;
   elseif a(1)==0, ini=1;
   elseif max(abs(roots(a))) > r, ini=1;
   else
      if a(1)~=1, a=a/a(1); end
      A = freqz(a,1,om); b = lslevin(M+1,om,A.*D,W./(abs(A).^2));
   end
else, ini=1;
end
if ini,     % compute FIR solution
   a = [1;zeros(N,1)]; b = lslevin(M+1,om,D,W);
end
x = [a(2:N+1);b]; delta = [];

% iterate (outer loop)
while 1,

   % compute complex error, Jacobian, and objective function value
   A = EM(:,1:N+1)*a; B = EM(:,1:M+1)*b; H = B./A;
   E = srW.*(D - H); l2error = E'*E; if N<1, break; end
   vec1 = srW./A; vec2 = -vec1.*H;
   J = [vec2(:,ones(1,N)).*EM(:,2:N+1),vec1(:,ones(1,M+1)).*EM(:,1:M+1)];

Could some kind soul perhaps provide an explanation or an alternative expression of the final two lines? I can't quite wrap my head around this being the computation of a matrix of partial derivatives as it stands.

Thanks

---

EDIT:

Having had a little play around with implementing a numerical approach to the Jacobian calculation, I think I can get that approach to work reasonably well. However, I still struggle to grasp the thesis implementation, which is presumably more direct/efficient.

My numerical approach, for reference:

% Compute Jacobian numerically.
PERTURB_STEP = 1e-4;
numericalJacobian = zeros(numel(H), nDenominator + 1 + mNumerator);

for iXTerm = 1 : nDenominator + 1 + mNumerator

    % Perturb target x term.
    x1 = x;
    x2 = x;
    x1(iXTerm) = x(iXTerm) - PERTURB_STEP;
    x2(iXTerm) = x(iXTerm) + PERTURB_STEP;

    % Separate a/b coefficient sets.
    a1 = [1; x1(1:nDenominator)];
    b1 = x1(nDenominator + 1 : nDenominator + 1 + mNumerator);

    a2 = [1; x2(1:nDenominator)];
    b2 = x2(nDenominator + 1 : nDenominator + 1 + mNumerator);

    % Get H terms.
    A1 = zPowMinusN * a1;
    B1 = zPowMinusM * b1;
    H1 = B1 ./ A1;

    A2 = zPowMinusN * a2;
    B2 = zPowMinusM * b2;
    H2 = B2 ./ A2;

    % Fill jacobian.
    numericalJacobian(:, iXTerm) = (H2 - H1) / (2 * PERTURB_STEP);

end

% Weight Jacobian.
numericalJacobian = numericalJacobian .* rootErrorWeight;

Answer 1

The Jacobian is not computed numerically but analytically and then just evaluated. The frequency response of the IIR filter is

$$H(e^{j\omega})=\frac{b_0+b_1e^{-j\omega}+\ldots+b_Me^{-jM\omega}}{1+a_1e^{-j\omega}+\ldots+a_Ne^{-jN\omega}}=\frac{B(e^{j\omega})}{A(e^{j\omega})}\tag{1}$$

Now you need the derivative with respect to the filter coefficients:

$$\frac{\partial H(e^{j\omega})}{\partial a_l}=-\frac{B(e^{j\omega})}{A^2(e^{j\omega})}e^{-jl\omega}=-\frac{H(e^{j\omega})}{A(e^{j\omega})}e^{-jl\omega}\tag{2}$$ and $$\frac{\partial H(e^{j\omega})}{\partial b_l}=\frac{e^{-jl\omega}}{A(e^{j\omega})}\tag{3}$$

The matrix $\boldsymbol{J}$ just contains all the derivatives $(2)$ and $(3)$ evaluated at all frequency points $\omega_k$. So the elements of $\boldsymbol{J}$ are given by

$$J_{kl}=\begin{cases}\displaystyle -\frac{H(e^{j\omega_k})}{A(e^{j\omega_k})}e^{-jl\omega_k},& 1\le l\le N\\\displaystyle \frac{e^{-jl\omega_k}}{A(e^{j\omega_k})},& N\lt l\le N+M+1\end{cases}\tag{4}$$

The only difference between the matrix J in the code and the matrix $\boldsymbol{J}$ in the text (and in $(4)$) is that in the code I also incorporated the (square root of the) weight matrix, so all matrix elements are multiplied by the square root of the weight at the corresponding frequency.