Approximation

To this point, we have discussed signal representations and models as basic tools for signal processing. In the following modules, we discuss the actual application of these tools to tasks such as approximation and compression, and we continue to discuss the geometric implications.

Linear approximation One common prototypical problem in signal processing is to find the best linear approximation to a signal x. By “best linear approximation,” we mean the best approximation to x from among a class of signals comprising a linear (or affine) subspace. This situation may arise, for example, when we have a noisy observation of a signal believed to obey a linear model. If we choose an ℓ2 error criterion, the solution to this optimization problem has a particularly strong geometric interpretation. To be more concrete, suppose S is a K-dimensional linear subspace of RN. (The case of an affine subspace follows similarly.) If we seek s * : = arg min s ∈ S s - x 2 , standard linear algebra results state that the minimizer is given by s * = A T A x , where A is a K×N matrix whose rows form an orthonormal basis for S. Geometrically, one can easily see that this solution corresponds to an orthogonal projection of x onto the subspace S (see (a)).

The linear approximation problem arises frequently in settings involving signal dictionaries. In some settings, such as the case of an oversampled bandlimited signal, certain coefficients in the vector α may be assumed to be fixed at zero. In the case where the dictionary Ψ forms an orthonormal basis, the linear approximation estimate of the unknown coefficients has a particularly simple form: rows of the matrix A in are obtained by selecting and transposing the columns of Ψ whose expansion coefficients are unknown, and consequently, the unknown coefficients can be estimated simply by taking the inner products of x against the appropriate columns of Ψ. For example, in choosing a fixed subset of the Fourier or wavelet dictionaries, one may rightfully choose the lowest frequency (coarsest scale) basis functions for the set S because, as discussed in Linear Models from Low-Dimensional Signal Models , the coefficients generally tend to decay at higher frequencies (finer scales). For smooth functions, this strategy is appropriate and effective; functions in Sobolev smoothness spaces are well-approximated using linear approximations from the Fourier or wavelet dictionaries . For piecewise smooth functions, however, even the wavelet-domain linear approximation strategy would miss out on significant coefficients at fine scales. Since the locations of such coefficients are unknown a priority, it is impossible to propose a linear wavelet-domain approximation scheme that could simultaneously capture all piecewise smooth signals. As an example, (a) shows the linear approximation of the Cameraman test image obtained by keeping only the lowest-frequency scaling and wavelet coefficients. No high-frequency information is available to clearly represent features such as edges.

Nonlinear approximation A related question often arises in settings involving signal dictionaries. Rather than finding the best approximation to a signal f using a fixed collection of K elements from the dictionary Ψ, one may often seek the best K-term representation to f among all possible expansions that use K terms from the dictionary. Compared to linear approximation, this type of nonlinear approximation , utilizes the ability of the dictionary to adapt: different elements may be important for representing different signals. The K-term nonlinear approximation problem corresponds to the optimization s K , p * : = arg min s ∈ Σ K ∥ s - f ∥ p . (For the sake of generality, we consider general Lp and ℓp norms in this section.) Due to the nonlinearity of the set ΣK for a given dictionary, solving this problem can be difficult. Supposing Ψ is an orthonormal basis and p=2, the solution to is easily obtained by thresholding: one simply computes the coefficients α and keeps the K largest (setting the remaining coefficients to zero). The approximation error is then given simply by the coefficients that are discarded: ∥ s K , 2 * - f ∥ 2 = ∑ k > K α ˜ k 2 1 / 2 . When Ψ is a redundant dictionary, however, the situation is much more complicated. We mention more on this below (see also (b)).

Measuring approximation quality One common measure for the quality of a dictionary Ψ in approximating a signal class is the fidelity of its K-term representations. Often one examines the asymptotic rate of decay of the K-term approximation error as K grows large. Defining σ K ( f ) p : = ∥ s K , p * - f ∥ p , for a given signal f we may consider the asymptotic decay of σK(f)p as K→∞. (We recall the dependence of and hence on the dictionary Ψ.) In many cases, the function σK(f)p will decay as K-r for some r, and when Ψ represents a harmonic dictionary, faster decay rates tend to correspond to smoother functions. Indeed, one can show that when Ψ is an orthonormal basis, then σK(f)2 will decay as K-r if and only if α˜k decays as k-r+1/2 .

Nonlinear approximation of piecewise smooth functions Let f∈CH be a 1-D function. Supposing the wavelet dictionary has more than H vanishing moments, then f can be well approximated using its K largest coefficients (most of which are at coarse scales). As K grows large, the nonlinear approximation error will decayWe use the notation f(α)≲g(α), or f(α)=O(g(α)), if there exists a constant C, possibly large but not dependent on the argument α, such that f(α)≤Cg(α). as σK(f)2≲K-H. Supposing that f is piecewise smooth, however, with a finite number of discontinuities, then (as discussed in Sparse (Nonlinear) Models from Low-Dimensional Signal Models) f will have a limited number of significant wavelet coefficients at fine scales. Because of the concentration of these significant coefficients within each scale, the nonlinear approximation rate will remain σK(f)2≲K-H as if there were no discontinuities present . Unfortunately, this resilience of wavelets to discontinuities does not extend to higher dimensions. Suppose, for example, that f is a CH smooth 2-D signal. Assuming the proper number of vanishing moments, a wavelet representation will achieve the optimal nonlinear approximation rate σK(f)2≲K-H/2 , . As in the 1-D case, this approximation rate is maintained when a finite number of point discontinuities are introduced into f. However, when f contains 1-D discontinuities (edges separating the smooth regions), the approximation rate will fall to σK(f)2≲K-1/2 . The problem actually arises due to the isotropic, dyadic supports of the wavelets; instead of O(1) significant wavelets at each scale, there are now O(2j) wavelets overlapping the discontinuity. We revisit this important issue in Compression.Despite the limited approximation capabilities for images with edges, nonlinear approximation in the wavelet domain typically offers a superior approximation to an image compared to linear approximation in the wavelet domain. As an example, (b) shows the nonlinear approximation of the Cameraman test image obtained by keeping the largest scaling and wavelet coefficients. In this case, a number of high-frequency coefficients are selected, which gives an improved ability to represent features such as edges. Better concise transforms, which capture the image information in even fewer coefficients, would offer further improvements in terms of nonlinear approximation quality.

Finding approximations As mentioned above, in the case where Ψ is an orthonormal basis and p=2, the solution to is easily obtained by thresholding: one simply computes the coefficients α and keeps the K largest (setting the remaining coefficients to zero). Thresholding can also be shown to be optimal for arbitrary ℓp norms in the special case where Ψ is the canonical basis. While the optimality of thresholding does not generalize to arbitrary norms and bases, thresholding can be shown to be a near-optimal approximation strategy for wavelet bases with arbitrary Lp norms . In the case where Ψ is a redundant dictionary, however, the expansion coefficients α are not unique, and the optimization problem can be much more difficult to solve. Indeed, supposing even that an exact K-term representation exists for f in the dictionary Ψ, finding that K-term approximation is NP-hard in general, requiring a combinatorial enumeration of the ZK possible sparse subspaces . This search can be recast as the optimization problem α ^ = arg min ∥ α ∥ 0 s.t. f = Ψ α . While solving is prohibitively complex, a variety of algorithms have been proposed as alternatives. One approach convexifies the optimization problem by replacing the ℓ0 fidelity criterion by an ℓ1 criterion α ^ = arg min ∥ α ∥ 1 s.t. f = Ψ α . This problem, known as Basis Pursuit , is significantly more approachable and can be solved with traditional linear programming techniques whose computational complexities are polynomial in Z. The ℓ1 criterion has the advantage of yielding a convex optimization problem while still encouraging sparse solutions due to the polytope geometry of the ℓ1 unit ball (see for example and ). Iterative greedy algorithms such as Matching Pursuit (MP) and Orthogonal Matching Pursuit (OMP) have also been suggested to find sparse representations α for a signal f. Both MP and OMP iteratively select the columns from Ψ that are most correlated with f, then subtract the contribution of each column, leaving a residual. OMP includes an additional step at each iteration where the residual is orthogonalized against the previously selected columns.

Manifold approximation We also consider the problem of finding the best manifold-based approximation to a signal (see (c)). Suppose that F={fθ:θ∈Θ} is a parametrized K-dimension manifold and that we are given a signal I that is believed to approximate fθ for an unknown θ∈Θ. From I we wish to recover an estimate of θ. Again, we may formulate this parameter estimation problem as an optimization, writing the objective function (here we concentrate solely on the L2 or ℓ2 case) D ( θ ) = f θ - I 2 2 and solving for θ * = arg min θ ∈ Θ D ( θ ) . We suppose that the minimum is uniquely defined. Standard nonlinear parameter estimation tells us that, if D is differentiable, we can use Newton's method to iteratively refine a sequence of guesses θ(0),θ(1),θ(2),⋯ to θ* and rapidly convergence to the true value. Supposing that F is a differentiable manifold, we would let J = [ ∂ D / ∂ θ 0 ∂ D / ∂ θ 1 ⋯ ∂ D / ∂ θ K - 1 ] T be the gradient of D, and let H be the K×K Hessian, Hij=∂2D∂θi∂θj. Assuming D is differentiable, Newton's method specifies the following update step: θ ( k + 1 ) ← θ ( k ) + [ H ( θ ( k ) ) ] - 1 J ( θ ( k ) ) . To relate this method to the structure of the manifold, we can actually express the gradient and Hessian in terms of signals, writing D ( θ ) = ∥ f θ - I ∥ 2 2 = ∫ ( f θ - I ) 2 d x = ∫ f θ 2 - 2 I f θ + I 2 d x . Differentiating with respect to component θi, we obtain ∂ D ∂ θ i = J i = ∂ ∂ θ i ∫ f θ 2 - 2 I f θ + I 2 d x = ∫ ∂ ∂ θ i ( f θ 2 ) - 2 I ∂ ∂ θ i f θ d x = ∫ 2 f θ τ θ i - 2 I τ θ i d x = 2 〈 f θ - I , τ θ i 〉 , where τθi=∂fθ∂θi is a tangent signal. Continuing, we examine the Hessian, ∂ 2 D ∂ θ i ∂ θ j = H i j = ∂ ∂ θ j ∂ D ∂ θ i = ∫ ∂ ∂ θ j 2 f θ τ θ i - 2 I τ θ i d x = ∫ 2 τ θ i τ θ j + 2 f θ τ θ i j - 2 I τ θ i j d x = 2 〈 τ θ i , τ θ j 〉 + 2 〈 f θ - I , τ θ i j 〉 , where τθij=∂2fθ∂θi∂θj denotes a second-derivative signal. Thus, we can interpret Newton's method geometrically as (essentially) a sequence of successive projections onto tangent spaces on the manifold. Again, the above discussion assumes the manifold to be differentiable. Many interesting parametric signal manifolds are in fact nowhere differentiable — the tangent spaces demanded by Newton's method do not exist. However, in we have identified a type of multiscale tangent structure to the manifold that permits a coarse-to-fine technique for parameter estimation.