Realizable algorithm for approximating Hilbert-Schmidt operators via Gabor Multipliers.

In this work, we consider new computational aspects to improve the approximation of Hilbert-Schmidt operators via generalized Gabor multipliers. One aspect is to consider the approximation of the symbol of an Hilbert-Schmidt operator as L2 projection in the spline-type space associated to a Gabor multiplier. This gives the possibility to employ a selection procedure of the analysis and synthesis function, interpreted as time-frequency lag; hence, with the related algorithm it is possible to handle both underspread and overspread operators. In the numerical section, we exploit the case of approximating overspread operators having compact and smooth spreading function and discontinuous time-varying systems. For the latter, the approximation of discontinuities in the symbol is not straightforward achievable in the generalized Gabor multipliers setting. For this reason, another aspect is to further process the symbol through a Hough transform, to detect discontinuities and to smooth them using a new class of approximants. This procedure creates a bridge between features detections techniques and harmonic analysis methods and in specific cases it almost doubles the accuracy of approximation.


Introduction
The problem of approximating Hilbert-Schmidt (HS) via Gabor multipliers was approached by several authors [4,8,11] due to its important applications in computational sciences and signal processing [1]. In [4], the authors have improved the approximation of HS operators using Generalized Gabor Multipliers (GGM). In [8], the authors present the connection between the approximation of the HS operators and GGM via multi-window spline-type ‡ The authors gratefully acknowledge the support of the Austrian Science Fund (FWF): project number P27516. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. spaces. The connection was presented only theoretically in that paper and no implementation was given. In this work, we have started from the theory proposed in [8] and we have designed a realizable algorithm that improves the overall numerical approximation of the HS operators in three ways: in speed by transferring the problems in splines-spaces domain, in complexity by lowering the number of atoms through a proper selection procedure, and geometrically via a feature detection technique i.e. the Hough transform. In a previous work [19], we have used the reformulation of Gabor system in spline-type spaces to produce a sparse-frequency representation of a signal. We use in this paper a generalized Gabor multipliers for the same purpose: selection of proper channels for the approximation of overspread or underspread operators. To do this, we will not use Euclidean spaces and we prefer an approach that uses locally compact(LC) groups as done in [14]. Also, we prefer a more classical approach [20] that does not involve modulation spaces [5] since we restrict our study to the class of HS operators. This general theory will lead to usual L 2 -projection once restricted to square integrable functions. The channels will not be orthogonal, but biorthogonal. For the class of HS operators we will obtain L 2 -projection, hence it can be seen as optimal design of biorthogonal frequency division multiplex (BFDM) [21] in a Gaborian dictionary.
In the Section 2 we outline a minimal necessary theory for introducing the novel concepts. In Section 3 we introduce GGM seen as elements of a spline-type (ST) space. In Subsection 3.1 we exploit the connection between GGM and ideals of the related ST spaces. This gives us a novel procedure for the selection of analysis and synthesis functions of a GGM. In Section 4 we will test the algorithm and analyse the outcome. The outcome is similar to the one developed in [7] but it handles the multi-window formulation of GGM. We show the robustness of the generators selection procedure for operators having smooth symbol and the inability of GGM of handling discontinuities. We overcome this problem by looking at the reminder as an image and using a common tool in image processing and pattern recognition: the Hough transform (HG) [15].

Theoretical brief
Given a LC group G the convolution for the space L 1 (G) can be defined as f * g ≔ ∫ G f(x)L x g(y)dx, where L y f(x) ≔ f(y −1 x) is the left shift. The morphisms from G into the torus form a LC group Ĝ called the dual group and each element, called character, can be represented as x ↦ 〈x, x〉 x ∈ G. The multiplication of a function with a character is called modulation, as in the case of the real groups ℝ d having as characters the Fourier basis 〈x, x〉 ≔ e 2iπx·x̂, x,x̂ ∈ ℝ d . The Fourier transform of a function in L 1 (G) is defined as x̂ ∈ Ĝ which composes with shift and character multiplication as (1) Onchis  (2) The "hat" notation is used in the mathematical literature for both the Fourier transform of function as well as for members of the dual group: ℱ(f)(x) or alternately f(x). Given σ ∈ S′ (ℝ G×Ĝ ), a tempered distribution, we can define the operator PDO σ : S(G) → S′(G) acting on the Schwartz class S(G) as (3) It is in fact an application of the Fourier transform to f, followed by multiplication with the mask σ, also called the Kohn-Nirenberg (K-N) symbol of the operator, and finally an application of the inverse Fourier transform. The K-N quantization rule establishes a correspondence between L 2 (G × Ĝ) and HS operator defined on L 2 (G) [22]. Another important element to analyse such an operator is its spreading function, which is the Fourier transform of its K-N symbol.The spreading function was mainly used in pseudodifferential operators theory [23,3] for the possibility of expressing an operator PDO σ as Our use of spreading function is instead closer to classification of operators. A PD (pseudodifferential) operator having spreading function well localized in 0 is called underspread, otherwise it is called overspread. The same distinction is made in the study of nonstationary random processes [17,16]. Overspread processes are related to high-lag time-frequency correlation; in Subsection 3.1 we will find the same lag behaviour analysing GGM through their spreading functions.

Gabor multipliers as Rihaczek spline-type systems
We want to produce an expansion of the symbol σ through well localized blocks in the phase space, as it happens for the approximation of a function through Gabor systems [10]. Given g ≔ {g i : i = 1,…, r} and γ ≔ {γ i : i = 1,…, r}, two sets of tempered distributions over G, and H a lattice in G × Ĝ, a Generalized Gabor multiplier (GGM for short) is an operator in the form (4) where m i is the mask and the left shift operator and modulation operator, usually denoted as L y and M x̂ [12], are joined in an unique time-frequency shift π(x,x) ≔ M x̂Ly . If the analysis and synthesis sets are composed by only one function the operator is called Gabor Multiplier (GM).
The rank one operator Q g,γ f ≔ 〈f, g〉 γ can be seen as a PD operator having the symbol called Rihaczek distribution of g and γ. The central property of this operator is that with L G×Ĝ being the shift that acts componentwisely [14], hence it is possible to avoid the inversion of the twisted convolution [13] that traditionally occurs by analysing the symplectic nature of the phase space [12]. By linearity of quantization, the approximation of a PD operator in the class of GGM generated by (g, γ, H) is a spline-type (ST) approximation of its symbol [6].

Definition 1.
Let G be a LC group, H be a closed and normal subgroup of G, be a Banach space, and be any set of functions or distributions in a translation invariant Banach space (ℬ, ∥ · ∥ ℬ ). We call spline-type system generated by Φ and H the collection of left shifts {L a ϕ i : a ∈ H, i ∈ I}, and spline-type space its closed span in ℬ.
In this paper we will consider only complex valued functions, i.e. .
In order to have a constructive realization of signal representation in a spline-type space we require the system to be a projection basis: given a Banach space B, a biorthogonal system in B × B′ is a family (ϕ i , ϕĩ) i∈I such that ϕĩ 2 (ϕ i l ) = δ i 1 ,i 2 . A biorthogonal system is a projection basis in X 0 ⊂ X if it is a basis for X 0 and (5) In formula (5) we can identify two operators: an analysis and a synthesis corresponding to the sets {ϕĩ} and {ϕ i }. The first is while the synthesis operator for {ϕ i } is U {ϕ} c(x) ≔ ∑ i∈I c i ϕ i , applied to a sequence c ≔ (c i ) i∈I .
Being a projection basis is equivalent for the set {ϕ i } to own a bounded synthesis [9]. Then a biorthogonal set exists and reproduction formula (5) holds.
In case of finitely generated ST spaces over L 2 (G) analysis and synthesis become and respectively, where c belongs usually to some (L p (H)) R space. The projection formula (5), can be rewritten for ST space as (6) Restricting our research to a discrete subgroup, we can express a characterization of projection basis, i.e. invertibility of the synthesis operator, as linear independence of the Fourier transforms of the generating set over shifts of the orthogonal of the lattice H. Onchis  The following statements are equivalent: ∄ξ ∈ Ĝ such that the functions on H ⊥ are linearly dependent This theorem is particularly appealing for many reasons:

(I)
The invertibility of the synthesis says that it is possible to identify S(Φ, H) with ideals of (L 1 (H),*).

(II)
In the proof of (iii) ⇒ (iii) an explicit algorithm BIORTHOGONAL for the construction of the biorthogonal system is stated.
(III) For symbols defined in L 2 (G×Ĝ) the projection (6) is the L 2 -projection [2]. We have best approximation of HS operators in GGM class .
By comparing (4) with the projection we get Proposition 2. Let g be compactly supported and γ having compactly supported Fourier transform. If its masks are of the form: It is important to notice that, since no symplectic computation is involved, the convolution in (7) can be quickly computed through the standard FFT.

GGM in the Fourier domain: selection procedure
One of the strength-weakness of a ST space is that the spectra of its generators defines most of its features. As pointed out in (I) we reproduce ideals of the convolution algebra (L 1 (H), *), which are in one-to-one correspondence with the ideals of the multiplication algebra (ℱL 1 (H), ·) through Fourier transform [20]. Informally, given a function f and a ST space generated by an atom g, if the spectra of f and g do not intersect, P g,H f = 0. Hence, to build a good ST space for the L 2 approximation of a function, it is enough to cover in a patch-like fashion its Fourier transform fulfilling requirement (iii) of Theorem 1. Straightforward question in our investigation is: What does each of this patch represent in the case of Gabor approximation? This quest will lead us in understanding the decomposition of a PD overspread operator by the particular GGM dictionary built through an analysis synthesis couple and some lagged versions. By only using properties (1) and (2) we obtain Straightforward drawback from the GGM construction is the inability of approximating symbols possessing important Fourier feature badly intersecting these patches. This will be shown in Section 4, when time varying filters will be discussed. We summarize the developed theory in the Algorithm 1.

Numerical tests
In numerical tests we consider HS operators defined over the class of discrete signals of finite length L as function over a periodic time domain ℤ L . As analysis and synthesis we have considered Gaussians or bump function as analysis and synthesis being the IFT of a bump function, in such a way that the related Rihaczek distribution is compactly supported, well localized in 0, and the SF is localized in 0 and quickly decaying. With M σ we will denote the obtained GGM approximation of the HS operator having symbol σ. We approach two types PD of operators: overspread operators and time-variant filters (TVF). The first shows the outcome of the dual shifts selection expressed in Subsection 3.1, while the latter exploit intrinsic limitation of the GGM approximation procedure. To build an overspread operator we have synthesized its SF σô ver as sum of tensor products 1D bumps. After fixing two Gaussian functions as analysis and synthesis

Algorithm 1 Selection of the generating set through Proposition 3 and computation of the GGM masks as in (7)
Precondition: Symbol σ P , Analysis window g, Synthesis windows γ, subgroup H, Shifts K̂ ⊂ Ĝ × G 17: end function window, we have selected a tolerance under which the SF of the GGM is considered null and selected the shifts over a rectangular lattice that intersect σô ver for at least of 1/5 of its support; we call this set Φ. The result is shown in Figure 1a. The SF of the GGM approximation is shown in Figure 1b. Computing the error in Frobenius norm, hence in HS norm and applying this GGM approximation M σ over , for a mean of 1000 random signal f we get Our interest in acoustics leads us in studying approximation of time-varying filters through GGM. In Figure 2a a simple example of non-stationary frequency filter is shown: its symbol σ filt takes only two values, 0 and 1, so it annihilates the frequency's content of the input signal starting from the high frequencies (the middle point of ). Since the symbol is characterized by a flat region, the most relevant harmonic is the 0th. For this reason we have considered a GM (i.e. mono-window case) generated by a Rihaczek distribution ϕ. We obtained E σ filt ≈ 0.0604, a quite small error predictable result in this 1-window case, but for a mean of 1000 randomly generated signal f we got E PDO σ filt ≈ 0.1224 with variance 5 · 10 −5 , which means a global relative error of the 12%. We need to analyse the reason for such a discrepacy between the two errors. As usual the answer is encapsulated in Fourier domain.
In Figure2b the magnitude of is shown in the Gaussian case. Low frequencies are not participating any longer, meaning that the GM approximates well the flat region of 1s, but lines running through 0 and having slope perpendicular to the jumps between 0 and 1 regions still persist. Using shifts selection does not help, and the error slightly increases. For this reason we decided to approximate these missed lines in 1D. We first employ a Hough transform to detect discontinuity in the symbol of the operator. From this, rotating by 90 we can have the lines l k on which the contribution in Fourier domain is relevant. This contribution can be then extracted and approximated by a favourable dictionary. We decide to use a reduced Gabor transform for fixed 1D Gaussian or bump atom, as done in [19]. After this correction we get which is a quite satisfactory outcome.

Conclusions
In this paper, we have introduced a realizable algorithm for HS operators approximation based on Gabor multipliers and multi-window spline-type spaces. We have used the spreading function of an operator in an operational way, getting to a novel interpretation of lag-behaviour of a channel through a Gabor dictionary. After the error analysis, we understood the limitation of this approach and come with a new use of the Hough transform for feature extraction over lines. 2a TVF and 2b SF reminder in GGM approximation Onchis and Zappalà Page 11