Levinson's algorithm and QR decomposition for complex least-squares FIR design

Question

I'm studying Mathias's thesis Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Response. In section 2.1.2 the complex LS approximation problem is defined by an overdetermined linear system $$ \mathbf{W}^{1/2}\mathbf{C}^H\mathbf{h}=\mathbf{W}^{1/2}\mathbf{d} \tag{2.8} $$ where $\mathbf{W}^{1/2}$ is a diagonal squared weighting matrix, $\mathbf{C}$ is a complex matrix which transforms the unknown impulse response vector $\mathbf{h}$ into frequency response, and $\mathbf{d}$ is the vector of desired response. This system can be represented by a normal equations $$ \mathbf{CWC}^H\mathbf{h} = \mathbf{CWd} \tag{2.9} $$

The unknown impulse response $\mathbf{h}$ can be solved by Eq. (2.9) as well as Eq. (2.8).

It is stated that solving Eq. (2.8) by QR decomposition is better from a numerical point of view because the condition number is squared in Eq. (2.9). Solving the normal equation, on the other hand, is better in terms of computational effort and memory requirement. This makes sense to me and I agree with it. So I try to compare these two different ways. Matt provides the Matlab code of the former and I write the second one. However the results don't support this view.

I use the example design in Matt's code: length 61 bandpass, band edges [.23,.3,.5,.57]*pi, weighting 1 in passband and 10 in stopbands, desired passband group delay 20 samples and another one in which passband group delay is 30 samples to make it linear phase while other requirements remains the same. The results are as follows.

Group delay = 20

Group delay = 30 (linear phase)

It is shown that Levinson's algorithm is better in both magnitude and phase approximation. The error norm of Levinson's in both cases is smaller than the one of QR decomposition. I also tried to solve the normal equation using QR decomposition by h3 = real(C * W * C') \ real(C * W * D); which results in the same filter coefficients as Eq. (2.8).

Another example is to approximates the magnitude response of a linear-phase audio EQ filter with coefficients

b = [1.04065742117985,-3.10743551019314,3.78294931146517,-2.20556080775822,0.510654247004686];
a = [1,-3.08533488056927,3.79819788644552,-2.22766143738208,0.536063093204188];

and the phase response is set to be linear. FIR length is 255 and the result shows in this case these two methods performs quite similarly but QR decomposition cannot give a strickly linear-phase filter.

So my question is

1. Why Levinson's algorithm outperforms QR decomposition in both terms of accuracy and efficiency?
2. What makes that the Levinson's algorithm able to design a strictly linear phase filter while QR decomposition cannot.

---

EDIT: Sorry I found something wrong with my code. The weighted error measure should be

    e = norm(Wsqrt * C' * h - Wsqrt * D);
    e2 = norm(Wsqrt * C' * h2 - Wsqrt * D);

After the bug fixed I found that in all cases the $\ell_2$ error of QR decomposition is indeed smaller than the error of Levinson's algorithm, which meets what Matt says in his thesis. So my questions should be then:

1. Is it a better choice to use QR decomposition instead of Levinson's algorithm for such problem since issue on computational effort and memory consumption is no longer a problem these days.
2. What makes that the Levinson's algorithm able to design a strictly linear phase filter while QR decomposition cannot.

---

Matlab code

N = 61;
groupdelay = 20;
om = pi * [linspace(0, .23, 230), linspace(.3, .5, 200), linspace(.57, 1, 430)];
D = [zeros(1, 230), exp(-1j * om(231:430) * groupdelay), zeros(1, 430)];
W = [10 * ones(1, 230), ones(1, 200), 10 * ones(1, 430)];
[h, h2, e, e2] = lslevin(N, om, D, W);

function [h, h2, e, e2] = lslevin(N, om, D, W)
    % h = lslevin(N,om,D,W)
    % Complex Least Squares FIR filter design using Levinson's algorithm
    %
    % h      filter impulse response
    % N      filter length
    % om     frequency grid (0 <= om <= pi)
    % D      complex desired frequency response on the grid om
    % W      positive weighting function on the grid om
    %
    % example: length 61 bandpass, band edges [.23,.3,.5,.57]*pi,
    % weighting 1 in passband and 10 in stopbands, desired passband
    % group delay 20 samples
    %
    % om=pi*[linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
    % D=[zeros(1,230),exp(-j*om(231:430)*20),zeros(1,430)];
    % W=[10*ones(1,230),ones(1,200),10*ones(1,430)];
    % h = lslevin(61,om,D,W);
    %
    % Author: Mathias C. Lang, Vienna University of Technology
    % 1998-07
    % [email protected]

    om = om(:); D = D(:); W = W(:); L = length(om);
    % DR = real(D); DI = imag(D);
    %% solve normal equation using Levinson's algorithm
    a = zeros(N, 1); b = a;

    % Set up vectors for quadratic objective function
    % (avoid building matrices)
    dvec = D; evec = ones(L, 1); e1 = exp(1j * om);

    for i = 1:N
        a(i) = W.' * real(evec);
        b(i) = W.' * real(dvec);
        evec = evec .* e1; dvec = dvec .* e1;
    end

    a = a / L; b = b / L;

    % Compute weighted l2 solution
    h = levin(a, b);
    
    %% solve original overdetermined linear system using QR decomposition
    n = (0:N-1).';   % building matrix
    C = cos(n*om.'); % real part of matrix C = exp(1j*n*om.')
    W = diag(W);
    Wsqrt = sqrt(W);
    h2 = (Wsqrt * C') \ real(Wsqrt * D);
    % h3 = real(C * W * C') \ real(C * W * D); % try to solve normal equation using QR decomposition

    %% weighted error measure
    e = norm(Wsqrt * C' * h);
    e2 = norm(Wsqrt * C' * h2);

end

function x = levin(a, b)
    % function x = levin(a,b)
    % solves system of complex linear equations toeplitz(a)*x=b
    % using Levinson's algorithm
    % a ... first row of positive definite Hermitian Toeplitz matrix
    % b ... right hand side vector
    %
    % Author: Mathias C. Lang, Vienna University of Technology, AUSTRIA
    % 1997-09
    % [email protected]

    a = a(:); b = b(:); n = length(a);
    t = 1; alpha = a(1); x = b(1) / a(1);

    for i = 1:n - 1
        k =- (a(i + 1:-1:2)' * t) / alpha;
        t = [t; 0] + k * flipud([conj(t); 0]);
        alpha = alpha * (1 - abs(k)^2);
        k = (b(i + 1) - a(i + 1:-1:2)' * x) / alpha;
        x = [x; 0] + k * flipud(conj(t));
    end

end

Answer 1

Nice work! However, I still believe that QR decomposition of the rectangular system matrix describing the overdetermined system $(2.8)$ is more numerically stable than the solution of the square system $(2.9)$. I'm not sure about the reason for your results, but I'm almost sure that there must be a bug somewhere. Both methods should give very accurate results for such a simple low order design problem.

I used Octave's implementation of the QR solver which is called when using the backslash operator to solve an overdetermined system (c.f. the corresponding doc-page).

I designed the filters in Octave using the following specs (which are the same as in your question):

w = pi * [linspace(0,.23,230),linspace(.3,.5,200),linspace(.57,1,430)];
D = [zeros(1,230),exp(-j*w(231:430)*20),zeros(1,430)];
W = [10*ones(1,230),ones(1,200),10*ones(1,430)];
w = w(:); D = D(:); W = W(:);
srW = sqrt(W);
N = 61;

The corresponding system $(2.8)$ was solved using the following code:

% solve overdetermined system using QR-decomposition
M = w * (0:N-1);
E = exp(-1i*M);
E = srW(:,ones(1,N)) .* E;
D = srW .* D;
h2 = [real(E);imag(E)] \ [real(D);imag(D)];

I solved the system $(2.9)$ using Levinson's algorithm as implemented by the function lslevin.m.

For the filter length $N=61$, both designs are identical, up to numerical precision.

Then I changed the filter length to $N=610$ (without changing the other specs), which gives a useless result because of the transition band which is much too wide for that filter length. But the idea was to create a system that is numerically problematic to solve.

Despite the absurd specs, both methods result in a relatively small approximation error (note that the transition band is a don't care band, which doesn't contribute to the error). But the $L2$-error obtained from solving $(2.9)$ using Levinson's algorithm was 2.74e-07, which is significantly larger than the approximation error 1.62e-10 resulting from solving $(2.8)$ using Octave's QR-solver.

This result seems to confirm my suspicion that in degenerate cases it is numerically more stable to directly solve the overdetermined system $(2.8)$ instead of the square system $(2.9)$. Luckily, in most practical cases I've encountered, with specs that are chosen carefully to match the desired filter length, the efficient solution of $(2.9)$ using Levinson's algorithm is very accurate.

---

Addressing the two questions in the edited part of your post:

1. I think that Levinson's algorithm gives sufficiently accurate results in almost all practical cases. It is efficient and it doesn't require storing a matrix due to the Toeplitz structure of the system of equations. I haven't yet come across a practical example for which Levinson's algorithm failed and where it was necessary to resort to QR decomposition to obtain a useful solution.

2. I haven't been able to reproduce the problem. If the desired group delay equals $(N-1)/2$, where $N$ is the number of taps, then the optimal solution has an exactly linear phase, and this solution can be computed either by Levinson's algorithm or by QR decomposition. In both cases, the resulting filter will have an exactly linear phase, up to numerical accuracy.