# Jacobian Computation in Least Squares IIR Filter Design

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
% 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;

• +1 for studying my thesis :) – Matt L. Jul 24 '19 at 9:52
• What is it that you don't understand, the derivation in the thesis, or its implementation in Matlab? – Matt L. Jul 24 '19 at 9:53
• Thanks for the quick response @MattL - and has to be good that it's the author himself! The MATLAB implementation I find most challenging, but I also didn't manage to complete a functional refactor from the derivation. I can see Eq 5.11, where at iteration k, you compute J from the gradient of H at x(k). I tried some central difference numerical approaches to this, but didn't get close to what your MATLAB implementation produces, unfortunately. – Nick M Jul 24 '19 at 10:03
• I'll have to take a look myself, it's not like it was yesterday ... :) – Matt L. Jul 24 '19 at 10:07
• No problem - thanks for the input :). I should also say that, though I evidently don't understand all of it, I'm really quite impressed by both the approach itself, and the performance of the filters I've produced while playing around with it. – Nick M Jul 24 '19 at 10:14

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.

• Awesome, thanks Matt. Equations 2 and 3 here in your response are where I'd gotten lost, so I'll have a go at doing the analytical partial derivative calculation, and avoid my clunky numerical version. – Nick M Jul 24 '19 at 19:57
• @NickM: You're welcome! Just use the chain rule and you'll arrive at the results (2) and (3). – Matt L. Jul 24 '19 at 21:03