A Douglas-Rachford construction of non-separable continuous compactly supported multidimensional wavelets
AA Douglas–Rachford construction of non-separablecontinuous compactly supportedmultidimensional wavelets
David Franklin ∗ Jeffrey A. Hogan ∗ Matthew K. Tam † June 8, 2020
Abstract
After re-casting the n -dimensional wavelet construction problem as a feasibilityproblem with constraints arising from the requirements of compact support, smooth-ness and orthogonality, the Douglas–Rachford algorithm is employed in the searchfor one- and two-dimensional wavelets. New one-dimensional wavelets are producedas well as genuinely non-separable two-dimensional wavelets in the case where thedilation on the plane is the standard D a f ( t ) = a − f ( t/a ) ( t ∈ R n , a > Keywords: wavelets; multiresolution analysis; optimisation; Douglas–Rachford algorithm; pro-jection algorithm; feasibility problem.
AMS subject classifications : 42C40, 42B99, 65T60, 47N10, 65K10, 65T60
Dedication
This paper is dedicated to the memory of Laureate Professor Jon Borwein, who first suggestedthis approach to the multidimensional wavelet construction problem. Jon was a friend andmentor to generations of mathematicians across the globe and has left an incomparable legacyof work spanning multiple disciplines. He was generous with his time and his ideas and was ahighly respected and well-loved faculty member at the University of Newcastle in Australia.
Continuous wavelet decompositions have been used in analysis since the 1930’s and in appliedmathematics since the 1980’s. They are implicit in the work of Calder´on on singular integrals[14] and explicit in the work of Grossman and Morlet on seismic exploration [26]. They may ∗ School of Mathematical and Physical Sciences, Mathematics Building, University of New-castle, University Drive, Callaghan NSW 2308,
Australia . ([email protected],jeff[email protected]) † School of Mathematics & Statistics, The University of Melbourne, Parkville VIC 3010,
Australia .([email protected]) a r X i v : . [ m a t h . C A ] J un e thought of as frame decompositions in which the index set associated with the frame is theupper half plane R = { ( x, a ) ∈ R : a > } , and the frame elements are generated from a single window function ψ by the action of dilations and translations. More precisely, given f ∈ L ( R ),we compute the frame coefficients W ψ f ( x, a ) by W ψ f ( x, a ) = (cid:104) f, ψ x,a (cid:105) = (cid:90) ∞−∞ f ( t ) 1 √ a ψ (cid:18) x − ta (cid:19) dt. (1)The mapping f (cid:55)→ W ψ f is known as the continuous wavelet transform (with respect to the wavelet ψ ). Given weak conditions on ψ , f may be recovered from the frame coefficients W ψ f ( x, a )( x ∈ R , a > ψ ) was developed [21], [30] – see also [25] for connectionswith the theory of singular integrals. Unfortunately, these constructions failed to generaliseto discrete data in such a way as to provide fast algorithms. On the other hand, Mallat [37]and Meyer [38] independently developed the concept of multiresolution analysis (MRA) whichenabled fast algorithms.Realisations of MRA’s require the construction of a scaling function ϕ ∈ L ( R ) with veryspecial properties. As a minimum, it is necessary that ϕ satisfies(i) { ϕ ( · − k ) } ∞ k = −∞ is an orthonormal collection in L ( R ).(ii) ϕ is self-similar in the sense that there exists a sequence { h k } ∞ k = −∞ ∈ (cid:96) ( Z ) such that12 ϕ (cid:16) x (cid:17) = ∞ (cid:88) k = −∞ h k ϕ ( x − k ) . (iii) (cid:82) ∞−∞ ϕ ( t ) dt = 1.Prototypical examples satisfying these conditions have long been known. The function ϕ H = χ [0 , , the characteristic function of [0 , Haar multiresolution analysis. Another example is ϕ S ( t ) = sin( πt ) / ( πt ) which is associated with the Shannon multiresolution analysis. Unfortunately, neither of these examples are satisfactory foruse in signal analysis and processing for reasons we outline below.When computing wavelet coefficients from (1), it is much preferred that the wavelet ψ becompactly supported, since this allows integration to be performed over a compact set. Infact, the shorter the support, the more efficiently this computation can be performed. Sincethe function ϕ S is not compactly supported (and, in fact, has very weak decay) it is thereforeunsuitable.The integral (1) represents time localised information about the signal f at scale a . With anapplication of the Parseval theorem for the Fourier transform, we have W ψ f ( x, a ) = (cid:90) ∞−∞ ˆ f ( ξ ) e − πixξ √ a ˆ ψ ( − aξ ) dξ (2)(where ˆ f and ˆ ψ are the Fourier transforms of f and ψ respectively). From (2) we see that thewavelet coefficients also give frequency localised information about ˆ f at the scale a − . For thisreason it is desirable that ˆ ψ also be compactly supported. Of course ψ and ˆ ψ cannot both becompactly supported, so we instead insist that ˆ ψ decay as fast as possible, or equivalently, that ψ be as smooth as possible. Hence, for the purpose of efficient numerics, we shall add the followingrequirements to the three conditions above: iv) ϕ is compactly supported.(v) ϕ is smooth.Note that the function ϕ H associated with the Haar multiresolution analysis fails condition (v),while the function ϕ S fails condition (iv). Without these properties, a multiresolution analysisfails to provide useful data and, in particular, without property (iv) a multiresolution analysiswill not provide fast algorithms for discrete data.Shortly after the publication of [37] and [38], Daubechies [18] used the MRA concept toconstruct a family of real-valued functions N ϕ which satisfy conditions (i)–(v) and for whichincreasing the support (indexed by the positive integer N ) gives improved smoothness. This ledto constructions of scaling functions ϕ with extra properties such as near-symmetry [19].Compactly supported wavelets with prescribed smoothness on R n can be easily generatedthrough tensor products of one-dimensional wavelets. However, such “separable” constructionssuffer from the preferential treatment of the directions associated with the coordinate axes, andproduce spurious artefacts in applications. Higher dimensional non-separable constructions haveproved elusive when one uses the obvious generalisation of the dilations suggested by the one-dimensional approach. On a more fundamental level, the one-dimensional constructions cannotbe easily transferred to higher dimensions as they involve techniques from complex analysissuch as spectral factorisations which are not available in multivariate complex analysis. Indeed,Kovaˇcevi´c and Vetterli [33] and Cohen and Daubechies [16] set out the theory of non-separablewavelets but did not explicitly construct any examples. Ayache [6] and Belogay and Wang [10]independently discovered methods of creating non-separable orthogonal wavelets in 1999. Theywere shortly followed by Lai and Roach [34], He and Lai [27] and Karoui [31, 32]. San Antolin andZalik [41] discovered a family of non-separable scaling functions and their associated framelets bymaking a change of variables in specific trigonometric polynomials. All of these methods generatenon-separable wavelets from one dimensional wavelets, typically by some kind of perturbation ormodulation. For a more detailed discussion of the methods used, we refer the reader to Lai [35]. Here we employ techniques from optimisation to construct new MRA-based one-dimensionalwavelets and new genuinely non-separable MRA-based multi-dimensional wavelets. We formulatethe design problem in terms of constraints on a matrix-valued function well-known to wavelettheorists, discretise the problem, and then numerically compute – through use of the
Douglas–Rachford algorithm – examples which simultaneously satisfy all of the constraints. This work isan extension of the PhD thesis of David Franklin [23]. A preliminary version of these resultsappears in [24].This paper is organised as follows. In Section 2 we review the basic axioms of a multiresolutionanalysis of L ( R n ) including details on how to encode properties of a scaling function ϕ into anassociated QMF m . These properties include the orthogonality of the integer shifts of ϕ , andthe compact support and regularity of ϕ . In Section 3, we consider the relevant constraints on m and the associated conjugate filters and express them in terms of constraints on a matrix-valued function U which has these filters as entries. We show that, in the case of compactlysupported scaling functions and wavelets, sampling can be used to discretise the constraints. InSection 4, the relevant background material in optimisation and the Douglas–Rachford algorithmfor solution of feasibility problems is introduced. This section provides a complete descriptionof the relevant Hilbert spaces, constraints and projections for the wavelet construction problem.Finally, Section 5 includes computational results of the application of the Douglas–Rachfordalgorithm to the one-dimensional and two-dimensional wavelet construction problems. .3 Notation We consider multi-indices α = ( α , α , . . . , α n ) ∈ Z n + , (i.e., each α i is a non-negative integer)and declare | α | = (cid:80) nj =1 α j . The partial order on multi-indices is defined by β ≤ α if and only if β j ≤ α j for 1 ≤ j ≤ n . By ∂ α we mean the differential operator ∂ α = (cid:18) ∂∂x (cid:19) α · · · (cid:18) ∂∂x n (cid:19) α n . The collection of N × N matrices with complex coefficients is denoted C N × N and the sub-collection of unitary matrices by U ( N ). The Frobenius norm of an N × N matrix A = ( a ij ) Ni,j =1 is given by (cid:107) A (cid:107) = (cid:16)(cid:80) Ni,j =1 | a ij | (cid:17) / .Given positive integers M and n , we define the set Q nM = { , , . . . , M − } n = { ( j , j , . . . , j n ) ∈ Z n ; 0 ≤ j i ≤ M − ≤ i ≤ n } . By ( C N × N ) Q nM we mean the collection of functions F : Q nM → C N × N . Elements of ( C N × N ) Q nM are known as matrix ensembles .The dot product of x , ξ ∈ R n is the real number (cid:104) x, ξ (cid:105) = (cid:80) nj =1 x j ξ j and we extend the dotproduct to z , ζ ∈ C n in the obvious way: (cid:104) z, ζ (cid:105) = (cid:80) nj =1 z j ζ j ∈ C .The Fourier transform ˆ f of f ∈ L ( R n ) is normalised by ˆ f ( ξ ) = (cid:82) R n f ( x ) e − πi (cid:104) x,ξ (cid:105) dx andextends unitarily to L ( R n ).A function f : R n → C is said to be Z n -periodic if f ( ξ + (cid:96) ) = f ( ξ ) for all ξ ∈ R n and (cid:96) ∈ Z n .The Lebesgue measure of a measurable subset E ⊂ R n is denoted | E | . The construction of a compactly supported smooth orthogonal scaling function–wavelet pair( ϕ, ψ ) on the line was first achieved by Daubechies in [18] with the help of the multiresolutionstructure introduced independently by Mallat [37] and Meyer [38]. The problem reduces to theconstruction of a periodic matrix-valued function U : R → C × satisfying certain restrictionsdesigned to force ϕ and ψ to have desirable properties for signal processing. The n -dimensionalwavelet construction problem may be reduced to the construction of a periodic matrix-valuedfunction U : R n → C n × n satisfying similarly motivated restrictions. The construction relies onthe notion of multiresolution analysis. In this section, we give an explanation of the multireso-lution structure and a discussion of the conditions we impose on the relevant filters to achievethese desirable properties. On L ( R n ) we have the unitary translation operators τ x ( x ∈ R n ) given by τ x f ( t ) = f ( t − x ). Let S be an n × n matrix with integer entries, all of whose eigenvalues have absolute value greater than1, and define an associated dilation operator D S on L ( R n ) by D S f ( t ) = (det( S )) − / f ( S − t ).There are of course many possibilities for the matrix S including (in two dimensions) the quincunxmatrix S = (cid:18) − (cid:19) . In this paper we consider only the matrices S = 2 I n (where I n is the n × n identity matrix) and in this case (with abusive notation) we write D = D I n . .2 Multiresolution analysis for L ( R n ) A multiresolution analysis ( { V j } ∞ j = ∞ , ϕ ) for L ( R n ) is a sequence of closed subspaces { V j } ∞ j = −∞ ⊂ L ( R n ) and a function ϕ ∈ V such that(i) V j ⊂ V j +1 for all j ∈ Z (ii) ∩ ∞ j = −∞ V j = { } and ∪ ∞ j = −∞ V j = L ( R n )(iii) f ∈ V j ⇐⇒ D − f ∈ V j +1 (iv) f ∈ V ⇐⇒ τ k f ∈ V ( k ∈ Z n )(v) { τ k ϕ } k ∈ Z n is an orthonormal basis for V . Orthonormality of the collection { τ k ϕ } k ∈ Z n is equivalent to the condition (cid:88) k ∈ Z n | ˆ ϕ ( ξ + k ) | = 1for almost every ξ . Given such a collection, we note that D ϕ ∈ V − ⊂ V and since { τ k ϕ } k ∈ Z n is an orthonormal basis for V , there exist constants { g k } k ∈ Z n ∈ (cid:96) ( Z n ) such that12 n ϕ (cid:16) x (cid:17) = (cid:88) k ∈ Z n g k ϕ ( x − k ) . (3)In fact, we have g k = 2 − n (cid:82) R n ϕ (cid:16) x (cid:17) ϕ ( x − k ) dx . Taking the Fourier transform of both sides of(3) gives ˆ ϕ (2 ξ ) = m ( ξ ) ˆ ϕ ( ξ ) (4)where m is the Z n -periodic Fourier series of { g k } , i.e., m ( ξ ) = (cid:80) k ∈ Z n g k e − πi (cid:104) k,ξ (cid:105) ( ξ ∈ R n ).Let V n be the vertices of the unit cube [0 , n in R n . Then | V n | = 2 n and if j ∈ { , , . . . , n − } has binary expansion j = (cid:80) n − k =0 a k k ( a k ∈ { , } ), we let v j = ( a , a , . . . , a n − ) ∈ V n .This provides a suitable enumeration of the elements of V n , i.e., V n = { v j } n − j =0 . Note that V = { , } ⊂ R and V = { v = (0 , , v = (1 , , v = (0 , , v = (1 , } ⊂ R . A necessary (but not sufficient) condition for the orthonormality of the collection { τ k ϕ } k ∈ Z n isthe quadrature mirror filter (QMF) condition n − (cid:88) j =0 | m ( ξ + v j / | = 1 (5)for almost every ξ .Since det(2 I n ) = 2 n , the index of the subgroup Z n / Z n is 2 n . Attached to each of the2 n − X ε of Z n / Z n (1 ≤ ε ≤ n −
1) is a subspace W ε and a waveletfunction ψ ε ∈ W ε such that V has the orthogonal decomposition V = V ⊕ W ⊕ W ⊕ · · · ⊕ W n − . (6) ith W εj = D j W ε we then have L ( R n ) = ⊕ ∞ j = −∞ ( ⊕ n − ε =1 W εj ) and the collection { j/ ψ ε (2 j x − k ) : j ∈ Z , k ∈ Z n , ≤ ε ≤ n − } forms an orthonormal basis for L ( R n ).Since D ψ ε ∈ W ε − ⊂ V , there are constants g εk such that12 n ψ ε (cid:16) x (cid:17) = (cid:88) k ∈ Z n g εk ϕ ( x − k ) (1 ≤ ε ≤ n − . (7)The Fourier transform of (7) may be written as (cid:99) ψ ε (2 ξ ) = m ε ( ξ ) ˆ ϕ ( ξ ), where m ε is the Z n -periodicFourier series of { g εk } , i.e., m ε ( ξ ) = (cid:80) k ∈ Z n g εk e − πi (cid:104) k,ξ (cid:105) ( ξ ∈ R n ). Given the orthonormality of { τ k ϕ } k ∈ Z n , the orthonormality of { τ k ψ ε } k ∈ Z n becomes equivalent to n − (cid:88) j =0 | m ε ( ξ + v j / | = 1 (8)for almost every ξ . Furthermore, the orthogonality of the decomposition (6) requires n − (cid:88) j =0 m ε ( ξ + v j / m η ( ξ + v j /
2) = δ εη (9)for almost every ξ . The requirement ∪ ∞ j = −∞ V j = L ( R n ) of a multiresolution analysis forces | ˆ ϕ (0) | = | (cid:82) R n ϕ ( t ) dt | =1. It is convenient to choose the phase of ϕ so that ˆ ϕ (0) = 1. Iterating equation (4) givesˆ ϕ ( ξ ) = m ( ξ/ m ( ξ/
4) ˆ ϕ ( ξ/
4) = · · · = J (cid:89) j =1 m ( ξ/ j ) ˆ ϕ ( ξ/ J ) . If m satisfies the QMF condition (5) and the infinite product (cid:81) ∞ j =1 m ( ξ/ j ) converges pointwisealmost everywhere, then its limit ˆ ϕ is square integrable and (cid:107) ϕ (cid:107) = 1 [19].It is relatively easy to see that if ϕ is supported on [0 , M − n ⊂ R n , then the coefficients h k = h k k ...k n in the dilation equation (3) are zero unless 0 ≤ k i ≤ M − ≤ i ≤ n ). Theconverse is trickier, and requires a higher-dimensional version of the Paley–Wiener theorem (seeTheorem 2.1 below).A function F : D ⊂ C n → C is holomorphic on D if for each z = ( z , z , . . . , z n ) ∈ D , thereis a polydisc P = { ( z , z , . . . , z n ) ∈ C n : | z − z | < r , . . . , | z n − z n | < r n } ⊂ D ( r , . . . , r n >
0) in which F may be represented by the absolutely convergent series F ( z ) = F ( z , z , . . . , z n ) = (cid:88) k ,k ,...,k n ≥ a k k ...k n ( z − z ) k ( z − z ) k · · · ( z n − z n ) k n . We say F is entire if it is holomorphic on D = C n . An entire function F : C n → C is of exponential type R > ε > A ε > | F ( z ) | ≤ A ε e π ( R + ε ) (cid:107) z (cid:107) here if z = ( z , z , . . . , z n ), (cid:107) z (cid:107) = (cid:80) nj =1 | z j | . The class of all functions of exponential type R > C n is denoted E n ( R ).Suppose F is the inverse Fourier transform of a function σ ∈ L ( R n ) which vanishes outside[ − R, R ] n = (cid:26) ξ ∈ R n : (cid:107) ξ (cid:107) ∞ = max ≤ j ≤ n | ξ j | ≤ R (cid:27) , i.e., F ( z ) = (cid:82) [ − R,R ] n σ ( ξ ) e πi (cid:104) z,ξ (cid:105) dξ where, if z = x + iy ∈ C n ( x, y ∈ R n ) and ξ ∈ R n , we have (cid:104) z, ξ (cid:105) = (cid:80) nj =1 z i ξ j = (cid:104) x, ξ (cid:105) + i (cid:104) y, ξ (cid:105) . Then F satisfies the pointwise bound | F ( z ) | ≤ (cid:90) [ − R,R ] n | σ ( ξ ) | e − π (cid:104) y,ξ (cid:105) dξ. However, if ξ ∈ [ − R, R ] n , then | ξ j | ≤ R for each 1 ≤ j ≤ n so that e − πy j ξ j ≤ e πR | y | and as aconsequence | F ( z ) | ≤ e πR (cid:107) y (cid:107) (cid:90) [ − R,R ] n | σ ( ξ ) | dξ ≤ n/ R n/ (cid:107) σ (cid:107) e πR (cid:107) z (cid:107) . Hence, F ∈ E n ( R ).The following multidimensional generalisation of the Paley–Wiener theorem is a special caseof a result given by Stein and Weiss [42] for more general support sets. Theorem 2.1 (Paley–Wiener theorem for cubes) . Suppose F ∈ L ( R n ) . Then F is the inverseFourier transform of a function vanishing outside the cube [ − R, R ] n if and only if F is therestriction to R n of a function in E n ( R ) . Given a positive integer N , we say Γ : C n → C is a trigonometric polynomial of degree N ifΓ( ζ ) = (cid:80) Nk ,k ,...,k n =0 a k k ...k n e − πi (cid:104) k,ζ (cid:105) ( ζ ∈ C n ) for some { a k k ...k n } Nk ,k ,...,k n =0 ⊂ C .The following result is a multi-dimensional version of Lemma 6.2.2 of [19]. Proposition 2.1.
Suppose Γ is a trigonometric polynomial of degree N on C n and Γ(0) = 1 .Let F ( ζ ) = ∞ (cid:89) j =1 Γ( ζ/ j ) ( ζ ∈ C n ) . Then F is the Fourier transform of a function f ∈ L ( R n ) supported on the cube [0 , N ] n .Proof. We prove the result in the case n = 2 only. Let Γ( ζ ) = (cid:80) Nk ,k =0 a k k e − πi (cid:104) k,ζ (cid:105) be as inthe statement of the proposition. Note that | e − πi (cid:104) k,ζ (cid:105) − | ≤ π (cid:107) k (cid:107) ∞ (cid:107) ζ (cid:107) so that | Γ( ζ ) | ≤ | Γ( ζ ) − | ≤ (cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) k ,k =0 a k k ( e − πi (cid:104) k,ζ (cid:105) − (cid:12)(cid:12)(cid:12)(cid:12) ≤ N (cid:88) k ,k =0 | a k ,k || e − πi (cid:104) k,ζ (cid:105) − |≤ N (cid:88) k ,k =0 π (cid:107) k (cid:107) ∞ (cid:107) ζ (cid:107) ≤ C (cid:107) ζ (cid:107) ≤ e C (cid:107) ζ (cid:107) (10) here C = 2 πN (cid:80) Nk ,k =0 | a k ,k | . However, if (cid:107) ζ (cid:107) ≤
1, from (10) we have (cid:12)(cid:12)(cid:12)(cid:12) (cid:81) ∞ j =1 Γ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:81) ∞ j =1 e C (cid:107) ζ (cid:107) / j = e C (cid:107) ζ (cid:107) ≤ e C , while if (cid:107) ζ (cid:107) > | e − πi (cid:104) k,ζ (cid:105) − | = | e − πi (cid:104) k,x (cid:105) e π (cid:104) k,y (cid:105) − e π (cid:104) k,y (cid:105) + e π (cid:104) k,y (cid:105) − |≤ e π (cid:104) k,y (cid:105) | e − πi (cid:104) k,x (cid:105) − | + | e π (cid:104) k,y (cid:105) − | ≤ y , y ≤
0. Here ζ = x + iy with x = ( x , x ), y = ( y , y ) ∈ R . Therefore, if (cid:107) ζ (cid:107) > y , y ≤
0, we have | Γ( ζ ) | ≤ | Γ( ζ ) − |≤ N (cid:88) k ,k =0 | a k ,k || e − πi (cid:104) k,ζ (cid:105) − | ≤ N (cid:88) k ,k =0 | a k ,k | = C. We choose an integer j ≥ j ≤ (cid:107) ζ (cid:107) < j +1 . Then (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 Γ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ j (cid:89) j =1 C ∞ (cid:89) j = j +1 e C (cid:107) ζ (cid:107) / j = C j exp( C (cid:107) ζ (cid:107) / j ) ≤ C (cid:107) ζ (cid:107) α with α = ln C ln 2 . We conclude that if y , y ≤ (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 Γ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ C max { , (cid:107) ζ (cid:107) α } . (11)Suppose now that y > y ≤
0. Then Γ( ζ ) = e − πiNζ ˜Γ( ζ ) with˜Γ( ζ ) = N (cid:88) k ,k =0 b k ,k e − πik ( − ζ ) e − πik ζ and b k ,k = a N − k ,k so that (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 Γ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) = | e − πiNζ | (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 | ˜Γ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ce πNy max { , (cid:107) ζ (cid:107) α } . (12)Similarly, if y ≤ y ≥ (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 Γ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ce πNy max { , (cid:107) ζ (cid:107) α } (13)and if y , y ≥ (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 Γ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ce πN ( y + y ) max { , (cid:107) ζ (cid:107) α } . (14)Let P ( ζ ) = e πiN ( ζ + ζ ) Γ( ζ ). Since | e πiN ( ζ + ζ ) | = e − πN ( y + y ) we have (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 P ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) = e − πN ( y + y ) ∞ (cid:89) j =1 | Γ( ζ/ j ) | . (15) pplying (15) to (11)–(14) gives (cid:12)(cid:12)(cid:12)(cid:12) ∞ (cid:89) j =1 P ( ζ/ j ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ce πN (cid:107) y (cid:107) max { , (cid:107) ζ (cid:107) α } ≤ C ε e π ( N + ε ) (cid:107) ζ (cid:107) for all ζ ∈ C n . By Theorem 2.1, the product (cid:81) ∞ j =1 P ( ζ/ j ) is the Fourier transform of a function σ ∈ L ( R n ) supported on the cube [ − N/ , N/ n . But F ( ζ ) = ∞ (cid:89) j =1 Γ( ζ/ j ) = e − πiN ( ζ + ζ ) ∞ (cid:89) j =1 P ( ζ/ j ) = e − πiN ( ζ + ζ ) ˆ σ ( ζ ) . If f ( x ) = σ ( x − ( N/ , N/ f is supported on [0 , N ] and ˆ f ( ζ ) = e − πiN ( ζ + ζ ) ˆ σ ( ζ ) = F ( ζ ). Since m satisfies (5), we have | m ( ξ ) | ≤ ξ . Since we require ˆ ϕ (0) = 1,equation (4) requires m (1) = 1. Because of the MRA condition (5), we also have m ( v j /
2) = 0for 1 ≤ j ≤ n −
1. The cross QMF condition (9) now gives m ε (0) = 0 for 1 ≤ ε ≤ n − m ( v j /
2) = δ j , m ε (0) = 0 (0 ≤ j ≤ n − , ≤ ε ≤ n − . (16) The following result is a consequence of [39, Chapter 3.7, Proposition 4].
Theorem 2.2.
Suppose ϕ is a compactly supported scaling function and { ψ ε } n − ε =1 is a collectionof wavelets associated with an MRA of L ( R n ) , all of which have bounded partial derivatives oforder less than or equal to d . Then(i) (cid:82) R n x α ψ ε ( x ) dx = 0 for | α | ≤ d and ≤ ε ≤ n − .(ii) The conjugate filters { m ε } n − ε =1 satisfy ∂ α m ε ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 = 0 for ≤ ε ≤ n − , | α | ≤ d . (17)Condition (ii) is not sufficient to ensure regularity of the wavelets { ψ ε } n − ε =1 . Nevertheless,this is the condition we impose in an attempt to enforce regularity, with the expectation that thelarger the value of d (i.e., the “flatter” the filters m ε (1 ≤ ε ≤ n −
1) at the origin) the higherthe regularity.Flatness of the coniugate filters at the origin as in (17) coupled with the cross QMF condition(9) gives the flatness of the quadrature filter m at the points { v j / } n − j =1 as the next result shows. Proposition 2.2.
Suppose { m ε } n − ε =0 are trigonometric polynomials satisfying (9) and (17) and m (0) = 1 . Then m satisfies ∂ α m ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ = v j / = 0 for ≤ j ≤ n − , | α | ≤ d . (18) .2.5 Non-separability It is a simple matter to construct smooth orthogonal compactly supported wavelets on R n througha tensor-product construction. If n = 2, we let ( ϕ , ψ ) and ( ϕ , ψ ) be smooth orthogonalcompactly supported one-dimensional scaling function-wavelet pairs and define a two-dimensionalscaling function Φ and three two-dimensional wavelets Ψ , Ψ , Ψ byΦ( ξ , ξ ) = ϕ ( ξ ) ϕ ( ξ ); (19)Ψ ( ξ , ξ ) = ϕ ( ξ ) ψ ( ξ ); Ψ ( ξ , ξ ) = ψ ( ξ ) ϕ ( ξ ); Ψ ( ξ , ξ ) = ψ ( ξ ) ψ ( ξ ) . Then the collection { Ψ ε } ε =1 generates a smooth orthogonal compactly supported wavelet basison R . Such systems, however, perform poorly in image processing applications, producingartefacts in the directions of the coordinate axes [33]. Here we seek non-separable wavelet basesin which neither the scaling function nor the wavelets can be decomposed as the tensor productof two functions of a single variable. Although non-separability is not imposed as a constraint,it is a simple matter to check whether scaling functions and wavelets generated by our methodsare separable.Suppose a two-dimensional scaling function Φ is supported on [0 , M ] and separable as in(19). Let m be the two-dimensional scaling filter associated with Φ and let m (1)0 , m (2)0 be theone-dimensional scaling filters associated with ϕ and ϕ respectively. Then m is separable: m ( ξ , ξ ) = m (1)0 ( ξ ) m (2)0 ( ξ )and since m (1)0 (0) = m (2)0 (0) = 1, we have m ( ξ , ξ ) = m ( ξ , m (0 , ξ ) . (20)Recalling that m ( ξ , ξ ) = (cid:80) M − k ,k =0 g k ,k e − πi ( k ξ + k ξ ) , (20) becomes M − (cid:88) k ,k =0 g k ,k e − πi ( k ξ + k ξ ) = M − (cid:88) k ,(cid:96) =0 g k ,(cid:96) e − πik ξ M − (cid:88) k ,n =0 g n,k e − πik ξ which is equivalent to g p,q = (cid:18) M − (cid:88) (cid:96) =0 g p,(cid:96) (cid:19)(cid:18) M − (cid:88) n =0 g n,q (cid:19) . (21)Let G be the M × M matrix with ( j, k )-th entry G j,k = g j,k (0 ≤ j, k ≤ M − G = ( G )(( G ) T ) T where = (1 , , . . . , T ∈ R M . As a measureof the separability of a two-dimensional scaling function ϕ , we compute its separability measure S ( ϕ ) = (cid:107) G − ( G )(( G ) T ) T (cid:107) where (cid:107) · (cid:107) is the Frobenius norm. Note that S ( ϕ ) = 0 if and only if ϕ is separable. We seekscaling functions with separability measure significantly larger than zero. Equations (5), (8) and (9) may be neatly organised as follows: the orthogonality of the collections { τ k ϕ } k ∈ Z n and { τ k ψ ε } k ∈ Z n (1 ≤ ε ≤ n −
1) and the orthogonality of the spaces they span requiresthat the matrix-valued Z n -periodic function U : R n → C n × n given by U ( ξ ) j,ε = m ε ( ξ + v j /
2) (0 ≤ j, ε ≤ n −
1) (22) s unitary for all ξ . When n = 1, U ( ξ ) is the 2 × U ( ξ ) = (cid:18) m ( ξ ) m ( ξ ) m ( ξ + ) m ( ξ + ) (cid:19) with ξ ∈ R , while when n = 2, U ( ξ ) is the 4 × m ( ξ , ξ ) m ( ξ , ξ ) m ( ξ , ξ ) m ( ξ , ξ ) m ( ξ + , ξ ) m ( ξ + , ξ ) m ( ξ + , ξ ) m ( ξ + , ξ ) m ( ξ , ξ + ) m ( ξ , ξ + ) m ( ξ , ξ + ) m ( ξ , ξ + ) m ( ξ + , ξ + ) m ( ξ + , ξ + ) m ( ξ + , ξ + ) m ( ξ + , ξ + ) with ξ = ( ξ , ξ ) ∈ R .The matrix-valued function U of (22) holds the key to our approach to wavelet constructionin one- and higher dimensions. In this section, we record the conditions on U which encodethe orthogonality, compact support and regularity conditions on the filters m ε of Section 2.2.1.Then we explore a sampling-based approach to discretisation of the problem through use of thediscrete Fourier transform. Finally, we use this discretisation to express the problem of waveletconstruction as a feasibility problem in which the constraint sets live in a finite-dimensionalHilbert space of matrix ensembles. In this section, the orthogonality, regularity and compact support conditions imposed on a scalingfunction and its associated wavelets are couched in terms of the matrix-valued function U of (22).It is clear from the form of (22) that there are strong relationships between the rows of U ( ξ ), andthese relationships – known here as consistency conditions – also must be accounted for whendesigning such matrices for wavelet construction. Let V n be as in Section 2.2.1. We endow V n with a group structure, thinking of it as ( Z ) n withcoordinate-wise addition modulo 2:( v j ⊕ v k ) (cid:96) = ( v j ) (cid:96) + ( v k ) (cid:96) (mod 2) . (23)Each j ∈ Y n = { , , . . . , n − } determines a permutation τ j of V n given by v τ j ( k ) = v j ⊕ v k and a permutation matrix σ j ∈ C n × n with ( k, (cid:96) )-th entry( σ j ) k(cid:96) = (cid:40) τ j ( k ) = (cid:96) v j ⊕ v k = v k ⊕ v j , σ j is symmetric. Further, since each m ε is Z n -periodic,[ σ j U ( ξ )] kε = n − (cid:88) (cid:96) =0 ( σ j ) k(cid:96) U ( ξ ) (cid:96)ε = (cid:88) { (cid:96) ; τ j ( k )= (cid:96) } U ( ξ ) lε = U ( ξ ) τ j ( k ) ,ε = m ε ( ξ + v τ j ( k ) / m ε ( ξ + ( v j ⊕ v k ) / m ε ( ξ + ( v j + v k ) /
2) = U ( ξ + v j / kε , rom which we conclude that U ( ξ + v j /
2) = σ j U ( ξ ) (25)for all ξ ∈ R n and all v j ∈ V n .Since addition in ( Z ) n is commutative, so too is the collection of matrices { σ j } j ∈ Y n . Proposition 3.1. If j ∈ Y n has binary representation j = (cid:80) n − k =0 a k k ( a k ∈ { , } ) then σ j decomposes as σ j = n − (cid:89) k =0 ( σ k ) a k . Proof.
Note that if j has binary representation as in the statement of the proposition, then v j = ( a , a , . . . , a n − ) = (cid:80) n − k =0 a k v k . We apply (25) repeatedly to find σ j U ( ξ ) = U ( ξ + v j / U (cid:18) ξ + n − (cid:88) k =0 a k v k / (cid:19) = ( σ n − ) a n − U (cid:18) ξ + n − (cid:88) k =0 a k v k / (cid:19) = n − (cid:89) k =0 ( σ k ) a k U ( ξ ) . (26)Since the matrices { σ j } j ∈ Y n commute, the product in (26) is independent of the order of thefactors, and the product is well-defined. Since U ( ξ ) is unitary, the result follows from (26). Corollary 3.1.
The consistency condition (25) holds for all j ∈ Y n and all ξ ∈ R n if and onlyif U ( ξ + v k /
2) = σ k U ( ξ ) for all ≤ k ≤ n − . The unitarity of the matrix U ( ξ ) of (22) is not sufficient to ensure the orthogonalities we require.In one dimension, Cohen’s condition [15] provides an easily checked sufficient condition. Thefollowing result (due to Bownik [13]) is a generalisation of the one-dimensional Cohen condition.
Theorem 3.1.
Suppose m ∈ C ∞ ( R n ) is Z n -periodic and is such that the infinite product ˆ ϕ ( ξ ) = (cid:81) ∞ j =1 m (2 − j ξ ) converges in L ( R n ) . Suppose also that there exists a compact set K ⊂ R n such that(i) K contains a neighbourhood of the origin;(ii) | K ∩ ( (cid:96) + K ) | = δ (cid:96) for all (cid:96) ∈ Z d ;(iii) m (2 − j ξ ) (cid:54) = 0 for all integers j > and all ξ ∈ K .Then { ϕ ( · − k ) } k ∈ Z n forms an orthonormal set. If m ∈ C D ( R n ) is Z n -periodic with D > n/ ,then the converse is true. Note that if m is a trigonometric polynomial, then it is infinitely differentiable. FromTheorem 3.1 we see that if m ε (0 ≤ ε ≤ n −
1) are trigonometric polynomials for which thematrix U ( ξ ) given by (22) is unitary, then the orthogonalities we require will be assured provided m has no zeroes on [ − / , / n . .1.3 Compact support In Section 2.2.2, we saw that ϕ being supported on [0 , M − n is equivalent to the Fourier series m being a trigonometric polynomial of the form m ( ξ ) = (cid:80) k ∈ Q nM h k e − πi (cid:104) k,ξ (cid:105) . The compactsupport of the wavelets ψ ε (1 ≤ ε ≤ n −
1) is equivalent to the Fourier series m ε (1 ≤ ε ≤ n − U = U ( ξ ) to also be a trigonometric polynomial: U ( ξ ) = (cid:88) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) (27)where for each k ∈ Q nM , A k is a constant 2 n × n matrix whose entries are the coefficients g εk . As in section 2.2.3, we note that the density of ∪ ∞ j = −∞ V j in L ( R n ) requires the conditions (16)on the Fourier series { m ε } n − ε . We define1 ⊗ U (2 n −
1) = (cid:26)(cid:18) T V (cid:19) : V ∈ U (2 n − (cid:27) . Here = (0 , , . . . , T ∈ C n . Since U ( ξ ) is unitary, the conditions (16) can be summarised as: U (0) ∈ ⊗ U (2 n − . We define C ⊗ C (2 n − × (2 n − = (cid:26)(cid:18) b T B (cid:19) : b ∈ C , B ∈ C (2 n − × (2 n − (cid:27) . To enable the regularity of the wavelets we construct, we impose condition (17) of Theorem 2.2on { m ε } n − ε =1 . As we saw in Corollary 2.2, this implies condition (18) on m . Together theseconditions may be written in terms of the matrix-valued function U of (22) as follows: ∂ α U ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 ∈ C ⊗ C (2 n − × (2 n − for 1 ≤ | α | ≤ d . (28)In summary, the problem of the construction of compactly supported orthogonal smoothscaling functions and wavelets on the line is equivalent to the following: Problem 3.1 (Scaling function/wavelet pairs in R n ) . Given an even integer M ≥ , we seekmatrices { A k } k ∈ Q nM ⊂ C n × n such that the trigonometric polynomial U : R n → C n × n givenby (27) satisfies the following three conditions:(i) U ( ξ ) is unitary for all ξ ∈ R n .(ii) U ( ξ + v j /
2) = σ j U ( ξ ) for all ξ ∈ R and ≤ j ≤ n − where σ j is as in (24).(iii) U (0) ∈ ⊗ U (2 n − .To allow for regularity of the associated scaling function/wavelet pairs we also impose(iv) ∂ α U ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 ∈ C ⊗ C (2 n − × (2 n − for ≤ | α | ≤ d . onditions (i)–(iv) do not guarantee the orthogonality of the integer shifts of the scalingfunction. Bownik’s sufficient condition for orthogonality may be written as follows:(v) U ( ξ ) (cid:54) = 0 for (cid:107) ξ (cid:107) ∞ = max ≤ j ≤ n | ξ j | ≤ / U ( ξ ) is the top left-hand entry of U ( ξ ). Our algorithms are designed to find examples ofsequences { A k } k ∈ Q nM for which the function U defined by (27) satisfies conditions (i)–(iv). Afterfinding such an example, we discard it if (v) is not satisfied. The assumption that the function U = U ( ξ ) is a trigonometric polynomial allows for discretisa-tion through sampling. We use this observation to recast conditions (i)-(iv) of Problem 3.1 intoconstraints on a finite number of coefficient matrices { A k } k ∈ Q nM .If B, C ∈ C N × N , we define the inner product (cid:104) B, C (cid:105) by (cid:104) B, C (cid:105) = tr( BC ∗ ) = (cid:80) Ni,j =1 b ij c ij .The norm arising from this inner product is the Frobenius norm (cid:107) · (cid:107) . Let L ([0 , n , C N × N )be the collection of measurable functions F : [0 , n → C N × N for which (cid:82) [0 , n (cid:107) F ( ξ ) (cid:107) dξ < ∞ .Given F, G ∈ L ([0 , n , C N × N ), we declare the inner product (cid:104) F, G (cid:105) to be (cid:104)
F, G (cid:105) = (cid:90) [0 , n (cid:104) F ( ξ ) , G ( ξ ) (cid:105) dξ. The sequence space (cid:96) ( Z n , C N × N ) is the collection of functions C : Z n → C N × N for which (cid:80) k ∈ Z n (cid:107) C k (cid:107) < ∞ . The inner product of B and C ∈ (cid:96) ( Z n , C N × N ) is given by (cid:104) B , C (cid:105) = (cid:80) k ∈ Z n (cid:104) B k , C k (cid:105) . The Fourier transform F : (cid:96) ( Z n , C N × N ) → L ([0 , n , C N × N ) given by F ( A )( ξ ) = (cid:80) k ∈ Z n A k e − πi (cid:104) k,ξ (cid:105) is a unitary mapping with inverse F − given by ( F − G ) k = (cid:82) [0 , n G ( ξ ) e πi (cid:104) k,ξ (cid:105) dξ whenever the integral converges. The space T NM,n of N × N matrix-valued trigonometric polyno-mials of degree less than M − T NM,n = (cid:26) P : R n → C N × N ; P ( ξ ) = (cid:88) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) for some { A k } k ∈ Q nM ⊂ C N × N (cid:27) and the finite sequence space X NM,n is given by X NM,n = { C ∈ (cid:96) ( Z n , C N × N ); C k = 0 if k / ∈ Q nM } . We note that T NM,n is a closed subspace of L ([0 , n , C N × N ) and X NM,n is a closed subspace of (cid:96) ( Z n , C N × N ). The Fourier transform may be restricted to X NM,n , and in doing so it becomes aunitary mapping of X NM,n onto T NM,n which we continue to denote F .The orthogonal projection P NM,n from L ([0 , n , C N × N ) onto T NM,n is given by P NM,n F ( ξ ) = (cid:90) [0 , n K NM ( ξ − η ) F ( ξ ) dη where K NM ( ξ ) = (cid:40) e − πi ( M − ξ + ··· + ξ n ) (cid:81) nj =1 sin( πMξ j )sin( πξ j ) if ξ (cid:54) = 0 M n if ξ = 0and the orthogonal projection R NM,n from (cid:96) ( Z n , C N × N ) onto X NM,n is given by( R NM,n ( C )) k = (cid:40) C k if k ∈ Q nM he sampling operator D NM,n : T NM,n → ( C N × N ) Q nM is given by ( D NM,n P ) j = P ( j/M ) ( j ∈ Q nM ) and there is an obvious isomorphism τ between ( C N × N ) Q nM and X NM,n , namely( τ A ) j = (cid:40) A j if j ∈ Q nM F M : ( C N × N ) Q nM → ( C N × N ) Q nM is given by( F M B ) j = (cid:88) k ∈ Q nM B k e − πi (cid:104) k,j (cid:105) /M , with inverse given by ( F − M A ) k = M − n (cid:80) j ∈ M n A j e πi (cid:104) j,k (cid:105) /M . Given U ∈ T NM,n , we form the matrix ensemble U ∈ ( C N × N ) Q nM by uniform sampling: the j -th entry of U is U j = U ( j/M )( j ∈ Q nM ), i.e., U = D NM,n U . Furthermore, if U ( ξ ) = (cid:80) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) , then U j = U ( j/M ) = (cid:88) k ∈ Q nM A k e − πi (cid:104) j,k (cid:105) /M ; A k = 1 M n (cid:88) j ∈ Q nM U j e πi (cid:104) j,k (cid:105) /M , (29)i.e., the ensembles U and A form a (finite) Fourier transform pair. For this reason, properties of U = U ( ξ ) may be encoded into its samples U j = U ( j/M ) ( j ∈ Q nM ) by way of its coefficients A k ( k ∈ Q nM ). Written in the “ensemble” notation, we denote the finite Fourier transform operationsof equation (29) as follows: U = F M A ; A = ( F M ) − U = ( F M ) − D NM U ; U = F A . These relationships are summarised in the following commuting diagram. A ∈ X NM,n U ∈ T NM,n U ∈ ( C N × N ) Q nM FF M D NM,n
Sampling and the discrete Fourier transform provide a means through which Problem 3.1 may bediscretised in the sense that the construction of a matrix-valued function U ( ξ ) ( ξ ∈ R n ) satisfyingthe conditions of Problem 3.1 may be replaced by the construction of finitely many matrices U j satisfying a compatible collection of conditions. The consistency condition (ii) of Problem 3.1 can be written in terms of the coefficient matrices { A k } k ∈ Q nM or the sampled matrices { U j = U ( j/M ) } j ∈ Q nM . Proposition 3.2.
Let P be the trigonometric polynomial P ( ξ ) = (cid:80) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) ( { A k } k ∈ Q nM ⊂ C n × n ) and P j = P ( j/M ) ( j ∈ Q nM ) . Then the following are equivalent:(i) P ( ξ + v (cid:96) /
2) = σ (cid:96) P ( ξ ) for all ξ ∈ R n and all (cid:96) ∈ { , , . . . , n − } (ii) σ (cid:96) A k = ( − k (cid:96) A k for all k = ( k , k , . . . , k n ) ∈ Q nM and all (cid:96) ∈ { , , . . . , n − } iii) P j + Mv (cid:96) / = σ (cid:96) P j for all j ∈ Q nM and all (cid:96) ∈ { , , . . . , n − } .Proof. Suppose P satisfies the consistency condition (i). Then (cid:88) k ∈ Q nM A k ( − k (cid:96) e − πi (cid:104) k,ξ (cid:105) = (cid:88) k ∈ Q nM A k e − πi (cid:104) k,ξ + v (cid:96) / (cid:105) = P ( ξ + v (cid:96) /
2) = σ (cid:96) P ( ξ ) = (cid:88) k ∈ Q nM σ (cid:96) A k e − πi (cid:104) k,ξ (cid:105) . Comparing coefficients in the sums on both sides of this equality gives σ (cid:96) A k = ( − k (cid:96) A k , hence(i) ⇒ (ii). A similar calculation gives the converse. Now suppose { A k } k ∈ Q nM satisfies (ii). Then P j + Mv (cid:96) / = (cid:88) k ∈ Q nM A k e − πi (cid:104) k,j + Mv (cid:96)/ (cid:105) /M = (cid:88) k ∈ Q nM A k ( − k (cid:96) e − πi (cid:104) j,k (cid:105) /M = σ (cid:96) (cid:88) k ∈ Q nM A k e − πi (cid:104) j,k (cid:105) /M = σ (cid:96) P j so that (ii) ⇒ (iii). The converse is proved similarly. The discretisation of the problem of constructing wavelet matrices U ( ξ ) in n dimensions relies onthe fact that the sampling operator D NM,n : T NM,n → X
NM,n is a multiple of a unitary operator. Aswe saw at the start of this section, the orthogonality of the collections { τ k ϕ } k ∈ Z n and { τ k ψ ε } k ∈ Z n (1 ≤ ε ≤ n −
1) and the orthogonality of the spaces they span requires that U ( ξ ) as given in(22) be unitary for all ξ . It is not sufficient to impose unitarity of the samples U j = U ( j/M ).To see this, consider the one-dimensional example U ( ξ ) = 14 [2( I + σ ) + (1 + i )( I − σ ) e − πiξ + (1 − i )( I − σ ) e − πiξ ] ( ξ ∈ R ) . Here I is the 2 × σ = (cid:18) (cid:19) . Since σ = I , U satisfies the consistencycondition U ( ξ + 1 /
2) = σU ( ξ ). Furthermore, U (0) = U (1 /
4) = I while U (1 /
2) = U (3 /
4) = σ ,all of which are unitary, yet U (1 /
8) = 12 ( I + σ ) which is not unitary since U (1 / ∗ U (1 /
8) =12 ( I + σ ). Proposition 3.3.
Let { A k } k ∈ Q nM ⊂ C n × n . The trigonometric polynomial U ( ξ ) = (cid:80) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) , ( ξ ∈ R n ) is unitary for all ξ ∈ R n if and only if U ( j/ (2 M )) is unitary for all j ∈ Q n M .Proof. Let J nM = { − M, . . . , , . . . , M − } n . If U is unitary for all ξ , then it is clearly unitaryat all ξ ∈ Q n M / (2 M ). Note that for all ξ ∈ R n , U ( ξ ) ∗ U ( ξ ) = (cid:88) m ∈ J nM B m e − πi (cid:104) m,ξ (cid:105) (30)with B m = (cid:80) k ∈ Q nM A ∗ k A m + k ( m ∈ J nM ). Suppose now that U is unitary at all points of Q n M / (2 M ), i.e., U ( j/ (2 M )) is unitary for all j ∈ Q n M . Then for all j ∈ Q n M we have I n = U ( j/ (2 M )) ∗ U ( j/ (2 M )) = (cid:88) m ∈ J nM B m e − πi (cid:104) m,j (cid:105) / (2 M ) . (31) y the orthonormality and completeness of the Fourier basis { e m } m ∈ Q n M (where e m ( j ) = e − πi (cid:104) m,j (cid:105) / (2 M ) ) in (cid:96) ( J nM , C ), we conclude from (31) that B m = δ m, I n and from (30) that U ( ξ ) ∗ U ( ξ ) = I n for all ξ ∈ R n .Note that Proposition 3.3 involves the sampled ensemble { U ( j/ (2 M )) } j ∈ Q n M rather than { U ( j/M ) } j ∈ Q nM . In Section 3.4, we’ll see that since Q n M / (2 M ) = ∪ n − k =0 ( Q nM + v k / /M, the ensemble { U ( j/ (2 M )) } j ∈ Q n M can be computed from U = { U ( j/M ) } j ∈ Q nM , so that unitarityof U ( ξ ) at all ξ can be achieved by the imposition of appropriate conditions on { U ( j/M ) } j ∈ Q nM . The regularity condition (iv) of Problem 3.1 can be written in terms of the coefficient matrices { A k } k ∈ Q nM or the sampled matrices { U j = U ( j/M ) } j ∈ Q nM . Proposition 3.4.
Let P be the trigonometric polynomial P ( ξ ) = (cid:80) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) ( { A k } k ∈ Q nM ⊂ C N × N ), P j = P ( j/M ) ( j ∈ Q nM ) and α = ( α , . . . , α N ) ∈ Z N + . Then the following are equivalent:(i) ∂ α P ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 ∈ C ⊗ C ( N − × ( N − (ii) (cid:80) k ∈ Q nM k α A k ∈ C ⊗ C ( N − × ( N − (iii) (cid:80) j ∈ Q nM c αj P j ∈ C ⊗ C ( N − × ( N − where c αj = (cid:80) k ∈ Q nM k α e πi (cid:104) j,k (cid:105) /M .Proof. We have A k = 1 M n (cid:80) j ∈ Q nM P j e πi (cid:104) j,k (cid:105) /M . Furthermore, ∂ α P ( ξ ) = ( − πi ) | α | (cid:88) k ∈ Q nM k α A k e − πi (cid:104) k,ξ (cid:105) so that (cid:88) k ∈ Q nM k α A k = 1 M n (cid:88) j ∈ Q nM P j (cid:88) k ∈ Q nM k α e πi (cid:104) j,k (cid:105) /M = 1 M n (cid:88) j ∈ Q nM c αj P j with c αj as in the statement of the proposition.The wavelet construction problem has now been recast as follows: Problem 3.2.
Given an even integer M ≥ , we seek a matrix ensemble U ∈ ( C n × n ) Q nM suchthat(i) the matrix ensembles U (cid:96) = (cid:26) U (cid:18) jM + v (cid:96) M (cid:19)(cid:27) j ∈ Q nM (0 ≤ (cid:96) ≤ n − are unitary;(ii) U j + Mv (cid:96) / = σ (cid:96) U j for all j ∈ Q nM/ ;(iii) U ∈ ⊗ U (2 n − .To allow for regularity of the associated scaling function and wavelets, we also impose(iv) (cid:80) j ∈ Q nM c αj U j ∈ C ⊗ C (2 n − × (2 n − for ≤ | α | ≤ d where c αj = (cid:80) k ∈ Q nM k α e πi (cid:104) j,k (cid:105) /M . .4 Wavelet feasibility problem Let M be even, d a non-negative integer, and σ j be the permutation matrix of (24). We define( C n × n ) Q nM σ = { U ∈ ( C n × n ) Q nM : U j + Mv (cid:96) / = σ (cid:96) U j , ( j ∈ Q ( M/ n , ≤ (cid:96) ≤ n − } (32)to be the collection of σ -consistent ensembles. It is straightforward to verify that ( C n × n ) Q nM is a vector space over C under the usual componentwise operations, and ( C n × n ) Q nM σ is a vectorsubspace. Moreover, by Proposition 3.2 we have F − M ( C n × n ) Q nM σ = { A ∈ ( C n × n ) Q nM ; ( − k (cid:96) A k = σ (cid:96) A k , k ∈ Q nM } . (33)We note that in the Fourier-side description (33) of ( C n × n ) Q nM σ , the condition applies individ-ually to each of the matrices A k of the ensemble A rather than on certain pairs of matrices asin (32). Lemma 3.1.
Let U be the trigonometric polynomial U ( ξ ) = (cid:80) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) ( ξ ∈ R n ) with { A k } k ∈ Q nM ⊂ C n × n and U be the matrix ensemble with j -th term U j = U ( j/M ) ( j ∈ Q nM ) .Then U (( j + v (cid:96) / /M ) = ( F M χ (cid:96) F − M U ) j where ( χ (cid:96) A ) k = e − πi (cid:104) k,v (cid:96) (cid:105) /M A k ( k ∈ Q nM ) .Proof. Observe that U (( j + v (cid:96) / /M ) = (cid:88) k ∈ Q nM A k e − πi (cid:104) k,j + v (cid:96) / (cid:105) /M = (cid:88) k ∈ Q nM e − πi (cid:104) k,v (cid:96) (cid:105) /M ( F − M U ) k e − πi (cid:104) j,k (cid:105) /M = (cid:88) k ∈ Q nM ( χ (cid:96) F − M U ) k e − πi (cid:104) k,j (cid:105) /M = ( F M χ (cid:96) F − M U ) j . Lemma 3.1 allows us to rephrase Proposition 3.3 as follows:
Proposition 3.5.
Let U ( ξ ) = (cid:80) k ∈ Q nM A k e − πi (cid:104) k,ξ (cid:105) ( A k ∈ C n × n ) be a trigonometric polyno-mial, U ∈ ( C n × n ) Q nM be the matrix ensemble with j -th entry U j = U ( j/M ) and U ( (cid:96) ) = F M χ (cid:96) F − M U ∈ ( C n × n ) Q nM (0 ≤ (cid:96) ≤ n − . Then U ( ξ ) is unitary for all ξ ∈ R n if and only if the matrix ensembles { U ( (cid:96) ) } n − (cid:96) =0 are allunitary. It’s important to note that if U ∈ ( C n × n ) Q nM σ and { U j } j ∈ Q nM/ are unitary, then all entriesof U are unitary since for j ∈ Q nM/ , we have( U j + Mv (cid:96) / ) ∗ U j + Mv (cid:96) / = ( σ (cid:96) U j ) ∗ ( σ (cid:96) U j ) = U ∗ j σ ∗ (cid:96) σ (cid:96) U j = U ∗ j U j = I n . Therefore, when imposing unitarity on entries of an ensemble U ∈ ( C n × n ) Q nM σ , it is enough toimpose unitarity on the sub-ensemble { U j } j ∈ Q nM/ . .4.2 Consistency Given an ensemble V ∈ ( C n × n ) Q nM , we extend it to a periodic mapping V : Z n → ( C n × n ) Q nM by declaring V j + Mp = V j ( j ∈ Q nM , p ∈ Z n ). We consider translation operators τ k ( k ∈ Z n )acting on ( C n × n ) Q nM by ( τ k V ) j = V j + k ( j, k ∈ Z n ) . (34)If r is an integer with 0 ≤ r ≤ n −
1, we define an operator T r on ( C n × n ) Q nM by( T r V ) j = 2 M M − (cid:88) m r =0 V ( j ,...,j r − ,m r ,j r +1 ,...,j n ) − e πi ( j r − m r − / /M . (35)If (cid:96) ∈ Y n has binary expansion (cid:96) = (cid:80) n − r =0 a r r ( a r ∈ { , } ) then we define T (cid:96) = n − (cid:89) r =0 ( T r ) a r . (36)Since the operators { T r } n − r =0 commute, the product in (36) is well-defined. Lemma 3.2.
As operators acting on periodisations of ensembles in ( C n × n ) Q nM , τ m and T (cid:96) ( m ∈ Z n , (cid:96) ∈ Y n ) commute, i.e., τ m T (cid:96) = T (cid:96) τ m .Proof. If 0 ≤ r (cid:54) = p ≤ n − s ∈ Z , then( τ sv p T r V ) j = 2 M M − (cid:88) m r =0 V ( j ,...,j r − ,m r ,j r +1 ,...,j p + s,...,j n ) − e − πi ( j r − m r +1 / /M = ( T r τ sv p V ) j while if 0 ≤ r = p ≤ n − s ∈ Z ,( τ sv r T r V ) j = 2 M M − (cid:88) m r =0 V ( j ,...,j r − ,m r ,j r +1 ,...,j n ) − e − πi ( j r + s − m r +1 / /M = 2 M M − (cid:88) m r =0 V ( j ,...,j r − ,m r + s,j r +1 ,...,j n ) − e − πi ( j r − m r +1 / /M = ( T r τ sv r V ) j . We conclude that τ sv p T r = T r τ sv p for all 0 ≤ p, r ≤ n −
1. Hence, if (cid:96) ∈ Y n and T (cid:96) is definedas in (36), we have τ sv p T (cid:96) = τ sv p n − (cid:89) r =0 ( T r ) a r = n − (cid:89) r =0 ( T r ) a r τ sv p = T (cid:96) τ sv p . Finally, if m = (cid:80) n − p =0 s p v p ∈ Z n then τ m T (cid:96) = n − (cid:89) p =0 τ s p v p T (cid:96) = T (cid:96) n − (cid:89) p =0 τ s p v p = T (cid:96) τ m . Proposition 3.6.
Suppose U ∈ ( C n × n ) Q nM σ and U ( (cid:96) ) = F M χ (cid:96) F − M U for some ≤ (cid:96) ≤ n − .Then U ( (cid:96) ) satisfies the consistency condition, i.e., U ( (cid:96) ) ∈ ( C n × n ) Q nM σ . roof. If U ( (cid:96) ) j is the j -th component of U ( (cid:96) ) , then with j = ( j , . . . , j n ) ∈ Q nM , U j = U ( j ,...,j n ) and 0 ≤ r ≤ n − U (2 r ) j = ( F − M χ r F M U ) j = 1 M n (cid:88) k ∈ Q nM (cid:88) m ∈ Q nM U m e − πi (cid:104) m,k (cid:105) /M e πi (cid:104) k,j (cid:105) /M e − πik r /M = 1 M n (cid:88) m ∈ Q nM U m M − (cid:88) k r =0 e πik r ( j r − m r − / /M n (cid:89) p =1 p (cid:54) = r M − (cid:88) k p =0 e πik p ( j p − m p ) /M = 1 M (cid:88) m ∈ Q nM U m δ j − m · · · δ j r − − m r − δ j r +1 − m r +1 . . . δ j n − m n − e πi ( j r − m r − / /M −
1= 2 M M − (cid:88) m r =0 U ( j ,...,j r − ,m r ,j r +1 ,...,j n ) − e πi ( j r − m r − / /M = ( T r U ) j , (37)i.e., U (2 r ) = T r U with T r the convolution operator defined in (35). If (cid:96) = (cid:80) n − r =0 a r r ( a r ∈{ , } ), then T (cid:96) is defined as in (36) and by (37) we have U ( (cid:96) ) = F M χ (cid:96) F − M U = F M n − (cid:89) r =0 ( χ r ) a r F − M U = n − (cid:89) r =0 ( F M χ r F − M ) a r U = n − (cid:89) r =0 ( T r ) a r U = T (cid:96) U . Hence, if U ∈ ( C n × n ) Q nM σ is a consistent ensemble, (cid:96) ∈ Y n and 0 ≤ p ≤ n −
1, Lemma 3.2 gives U ( (cid:96) ) j + Mv p / = ( τ Mv p / T (cid:96) U ) j = ( T (cid:96) τ Mv p / U ) j = ( T (cid:96) σ p U ) j = σ p ( T (cid:96) U ) j = σ p U ( (cid:96) ) j , i.e., U ( (cid:96) ) is a consistent ensemble.We close this section with the discretised version of the wavelet construction problem: Problem 3.3.
Given an even integer M ≥ , we seek a matrix ensemble U ∈ ( C n × n ) Q nM σ suchthat(i) U ( (cid:96) ) = {F M χ (cid:96) F − M U } n − (cid:96) =1 are unitary ensembles;(ii) U ∈ ⊗ U (2 n − .To allow for regularity of the associated scaling function and wavelets, we also impose(iii) (cid:80) j ∈ Q nM c αj U j ∈ C ⊗ C (2 n − × (2 n − for ≤ | α | ≤ d where c αj = (cid:80) k ∈ Q nM k α e πi (cid:104) j,k (cid:105) /M . In this section we give the background required to solve Problem 3.3 with techniques borrowedfrom optimisation. .1 Projection operators Let H be a finite-dimensional Hilbert space. Given a set S ⊆ H , its (metric) projector is theset-valued operator given by P S ( x ) := { s ∈ S : (cid:107) s − x (cid:107) ≤ d ( x, S ) } ( x ∈ H )where d ( x, S ) = inf s ∈ S (cid:107) x − s (cid:107) . It is straightforward to check that P S ( x ) (cid:54) = ∅ for all x ∈ H solong as S is nonempty and closed. In a common abuse of notation, we write P S ( x ) = p to mean P S ( x ) = { p } . Proposition 4.1 (Properties of projectors) . Let H be a finite dimensional Hilbert space.(a) Let C , C , . . . , C m ⊆ H be nonempty closed sets and define C := C × · · · × C m ⊆ H m . Then P C = P C × · · · × P C m . (b) Let L : H → H be an isometric isomorphism and C ⊆ H be a nonempty closed set. Then P L ( C ) = L ◦ P C ◦ L − . Proof. (a): Follows easily from the definition.(b): Let x ∈ H . First note that since L is an isometric isomorphism, we have d ( x, L ( C )) = d ( L − x, C ). On one hand, if p ∈ P L ( C ) ( x ), then L − p ∈ C and d ( L − x, C ) = d ( x, L ( C )) = (cid:107) x − p (cid:107) = (cid:107) L − x − L − c (cid:107) . This implies that L − p ∈ P C ( L − x ) or, equivalently, that p ∈ ( L ◦ P C ◦ L − )( x ). On the otherhand, if p ∈ ( L ◦ P C ◦ L − )( x ), then there exists c ∈ P C ( L − x ) such that p = Lc and d ( x, L ( C )) = d ( L − x, C ) = (cid:107) L − x − c (cid:107) = (cid:107) x − Lc (cid:107) = (cid:107) x − p (cid:107) , which implies that p ∈ P L ( C ) ( x ). This completes the proof.In what follows, the unit sphere is denoted S := { x ∈ H : (cid:107) x (cid:107) = 1 } . We recall that thesingular value decomposition (SVD) of a matrix A ∈ C N × N is of the form A = U Σ V ∗ where U, V ∈ U ( N ) and Σ ∈ C N × N is a diagonal matrix with the diagonal entries (the singular values of A ) being the eigenvalues of √ A ∗ A . Proposition 4.2 (Examples of projectors) . Let H , H (cid:48) be finite dimensional Hilbert spaces.(a) Let L : H → H (cid:48) be linear and denote C := { x ∈ H : Lx = 0 } . If LL ∗ is invertible, then P C ( x ) = x − L ∗ ( LL ∗ ) − ( Lx ) ∀ x ∈ H . (b) Let x ∈ H . Then P S ( x ) = (cid:40) x (cid:107) x (cid:107) x (cid:54) = 0 , S x = 0 . (c) Let X ∈ C N × N . Then P U ( N ) ( X ) = { U V ∗ : X = U Σ V ∗ is an SVD } . Proof. (a): See [7, Example 28.14(iii)].(b): Follows easily from the definitions.(c): See [29, Theorems 8.1 & 8.6].We note that if σ ∈ U ( N ) and X ∈ C N × N then P U ( N ) ( σX ) = σP U ( N ) ( X ) . (38) .2 Projection Algorithms and Feasibility Problems Given finitely many closed sets C , . . . , C m ⊆ H (a finite-dimensional Hilbert space) with nonemptyintersection, the corresponding feasibility problem isfind x ∈ m (cid:92) k =1 C k . (39) Projection algorithms are a family of iterative algorithms which can be used to solve (39) by ineach step utilising only projectors onto the individual sets (rather than the entire intersectionat once). The two most important examples of projection algorithms are the method of cyclicprojections [12] and the
Douglas–Rachford (DR) method [36, 9], as well as their variants [11, 3].In this work we employ the
Douglas–Rachford method which can be compactly described asthe following fixed point iteration: Given x ∈ H , choose any sequence ( x k ) satisfying x k +1 ∈ T ( x k ) where T := I + R C R D , (40)and R A := 2 P A − I denotes reflector with respect to a set A . Here we note that the sequence( x k ) is only required to satisfy the inclusion in (40) since, in general, the operator T : H → H is a point-to-set mapping.When applying a method based on (40), the sequence of interest ( i.e., the one that solves(39)) is not ( x k ) itself, but one of its projections onto the set D . For this reason, it is convenientto implement the Douglas–Rachford algorithm as outlined in Algorithm 1 and, in order to beconcrete, we state a general convergence result for the convex setting in Theorem 4.1. Algorithm 1:
Implementation of the Douglas–Rachford algorithm.
Input: x ∈ H ;Set k := 0 and choose any p ∈ P D ( x ); while stopping criteria not satisfied do Choose any point x k +1 satisfying x k +1 ∈ x k + P C (2 p k − x k ) − p k ;Choose any point p k +1 satisfying p k +1 ∈ P D ( x k +1 );Set k := k + 1; endReturn: p k ; Although Algorithm 1 applies to problem (39) with n = 2, the general problem (39) canalways be cast as a two set problem via the following product space formulation. Let C , D besubsets of H m given by C := C × C × · · · × C m , D := { ( x, x, . . . , x ) ∈ H m : x ∈ H} . hen the following equivalence holds: x ∈ m (cid:92) k =1 C k ⇐⇒ ( x, x, . . . , x ) ∈ C ∩ D. From here onwards, when speaking of applying the Douglas–Rachford algorithm to a feasibilityproblem, we will always mean its product space reformulation.
Theorem 4.1 (Behaviour of the DR algorithm [9, Theorem 3.13]) . Suppose
C, D ⊆ H are closedand convex with nonempty intersection. Let x ∈ H and set x k +1 = T ( x k ) for all k ∈ N . Thenthe sequence ( x k ) converges to a point x ∈ Fix T := { x : T x = x } and, moreover, P D ( x ) ∈ C ∩ D . In general, beyond the case of convex sets there is insufficient theory to justify applicationof projection methods. Indeed, most non-convex results in the literature rely on restrictive regu-larity notions from nonsmooth analysis and, even then, only yield local convergence guarantees[28, 40, 17]. Nevertheless, projection methods have been empirically observed to still performreasonably well in certain non-convex settings include matrix completion [2], graph colouring[4], combinatorial optimization [5, 1], road design [8], and constraint satisfaction [22]. This ex-perience suggests use of the Douglas–Rachford method in the setting outlined in the followingsection.
Although the matrices we work with have complex entries, for the purpose of algorithms is moreconvenient to work in a space over the real field. In this section, we provide the necessarybackground to justify this process. Before doing so, we first recall that the
Frobenius inner-product on C N × N , denoted (cid:104)· , ·(cid:105) F , is given by (cid:104) U, V (cid:105) F := Tr( U ∗ V ) = (cid:80) Ni,j =1 U ∗ ij V ij . The inducednorm is known as the Frobenius norm and is given by (cid:107) U (cid:107) = N (cid:88) i,j =1 | U ij | = n (cid:88) i,j =1 ( (cid:60) U ij ) + ( (cid:61) U ij ) , (41)where (cid:60) z and (cid:61) z denote the real and imaginary parts of a complex number z , respectively.Given a finite set A with | A | = m , we consider the collection ( C N × N ) A of matrix-valuedfunctions F : A → C N × N which, with abusive notation, we identify with H := ( C N × N ) m := C N × N × · · · × C N × N (cid:124) (cid:123)(cid:122) (cid:125) m factors . Depending on the inner-product and field, ( C N × N ) m may be viewed as a Hilbert space in twoways:(a) Over the field C , H can be equipped with the inner-product (cid:104)· , ·(cid:105) C given by (cid:104) U , V (cid:105) C := m (cid:88) j =1 (cid:104) U j , V j (cid:105) F . (42)(b) Over the field R , H can be equipped with the inner-product (cid:104)· , ·(cid:105) R given by (cid:104) U , V (cid:105) R := (cid:104)(cid:60) U , (cid:60) V (cid:105) F + (cid:104)(cid:61) U , (cid:61) V (cid:105) F = m (cid:88) j =1 (cid:104)(cid:60) U j , (cid:60) V j (cid:105) F + m (cid:88) j =1 (cid:104)(cid:61) U j , (cid:61) V j (cid:105) F . (43) ince we will only be concerned with the latter (real) inner-product, we will drop the subscript“ R ” whenever there is no ambiguity. Proposition 4.3.
The norms in both of the aforementioned spaces coincide.Proof.
Follows by combining (41), (42) and (43).
We concentrate now on the Hilbert space H = ( C n × n ) Q nM σ and observe that the discretisedwavelet construction Problem 3.3 is equivalent to the following: Problem 4.1.
Given an integer M ≥ , find a matrix ensemble U = { U j } j ∈ Q nM ∈ ∩ n − (cid:96) =0 C ( (cid:96) )1 ∩ C ⊂ H where the constraint sets are defined as C (0)1 := { U ∈ H : U j ∈ U (2 n ) ( j ∈ Q nM/ \ { } ) , U ∈ ⊗ U (2 n − } C ( (cid:96) )1 := { U ∈ H : ( F − M χ (cid:96) F M U ) j ∈ U (2 n ) , ( j ∈ Q nM/ , ≤ (cid:96) ≤ n − } C := (cid:26) U ∈ H : (cid:88) k ∈ Q nM c αk U k ∈ C ⊗ C (2 n − × (2 n − for ≤ | α | ≤ d (cid:27) where c αk = (cid:80) j ∈ Q nM j α e πi (cid:104) j,k (cid:105) /M . C ( (cid:96) )1 Recall that by Proposition 3.5, unitarity of the trigonometric polynomial U ( ξ ) at all ξ is equiv-alent to the unitarity of the ensembles { U ( (cid:96) ) } n − (cid:96) =0 where U = { U j = U ( j/M ) } j ∈ Q nM and U ( (cid:96) ) = F M χ (cid:96) F − U . Completeness requires m (0) = 1 or equivalently, ( U ) = 1.Let A = (cid:18) a v T w B (cid:19) ∈ C n × n with a ∈ C , v , w ∈ C n − and B ∈ C (2 n − × (2 n − . Then theprojection P ⊗U (2 n − ( A ) of A onto 1 ⊗ U (2 n −
1) is given by P ⊗U (2 n − ( A ) = (cid:18) T P U (2 n − ( B ) (cid:19) , so the projection P (0) C from H onto C is given by( P (0) C U ) j = (cid:40) σ (cid:96) P ⊗U (2 n − ( U ) if j = M v (cid:96) / P U (2 n ) U j if j (cid:54) = M v (cid:96) / P U (2 n ) is the projection of C n × n onto U (2 n ) given in Proposition 4.2.We recall the translation operators τ k of equation (34) and define modulation operators µ k ( k ∈ Z n ) on ( C n × n ) Q nM defined by ( µ k U ) j = e πi (cid:104) j,k (cid:105) /M U j . We then have the intertwiningrelations F M τ k = µ k F M ; F M µ k = τ − k F M ( k ∈ Z n ) (44) nd similarly, F − M τ k = µ − k F − M and F − M µ k = τ k F − M . The relationship between the modulationoperators µ k and the operator χ (cid:96) of Lemma 3.1 is given by µ k = (cid:81) n(cid:96) =1 ( χ (cid:96) ) − k (cid:96) from which weimmediately see that µ k χ (cid:96) = χ (cid:96) µ k . (45)Let S (cid:96) = F M χ (cid:96) F − M (1 ≤ (cid:96) ≤ n ). Then (44) and (45) give S (cid:96) τ k = F M χ (cid:96) F − M τ k = F M χ (cid:96) µ k F − M = F M µ k χ (cid:96) F − M = τ k F M χ (cid:96) F − M = τ k S (cid:96) . (46)Let P U (2 n ) QnM be the projection of ( C n × n ) Q nM onto U (2 n ) Q nM = { U ∈ ( C n × n ) Q nM : U j ∈ U (2 n ) for all j ∈ Q nM } given by ( P U (2 n ) QnM U ) j = P U (2 n ) U j ( j ∈ Q nM ) . For 1 ≤ (cid:96) ≤ n −
1, consider the operator P C ( (cid:96) )1 : ( C n × n ) Q nM → ( C n × n ) Q nM given by P C ( (cid:96) )1 U = S − (cid:96) P U (2 n ) QnM S (cid:96) . We aim to show that P C ( (cid:96) )1 is the projection of H onto C ( (cid:96) )1 .Given A ∈ ( C n × n ) Q nM and X ∈ C n × n , we define X A ∈ ( C n × n ) Q nM by ( X A ) j = XA j ( j ∈ Q nM ). Note that if σ ∈ U (2 n ) and A ∈ ( C n × n ) Q nM then an application of (38) gives[ P U (2 n ) QnM ( σ A )] j = P U (2 n ) ( σA j ) = σP U (2 n ) ( A j ) = σ ( P U (2 n ) QnM A ) j = [ σ P U (2 n ) QnM A ] j so that P U (2 n ) QnM ( σ A ) = σ P U (2 n ) QnM A . Proposition 4.4.
Suppose U ∈ H is an ensemble satisfying the consistency condition, i.e., τ − Mv k / U = σ k U (0 ≤ k ≤ n − . Then for ≤ (cid:96) ≤ n − , P C ( (cid:96) )1 U ∈ H , i.e., P C ( (cid:96) )1 preserves H .Proof. Since τ n P U (2 n ) QnM = P U (2 n ) QnM τ n , we apply (46) and the consistency condition to find τ − Mv k / P C ( (cid:96) )1 U = τ − Mv k / S − (cid:96) P U (2 n ) QnM S (cid:96) U = S − (cid:96) τ − Mv k / P U (2 n ) QnM S (cid:96) U = S − (cid:96) P U (2 n ) QnM τ − Mv k / S (cid:96) U = S − (cid:96) P U (2 n ) QnM S (cid:96) τ − Mv k / U = S − (cid:96) P U (2 n ) QnM S (cid:96) σ k U = σ k S − (cid:96) P U (2 n ) QnM S (cid:96) U = σ k P C ( (cid:96) )1 U which completes the proof. C An ensemble U ∈ ∩ n − (cid:96) =1 C ( (cid:96) )1 may be interpreted as samples of a trigonometric polynomial U : R n → ( C n × n ) Q nM σ . In fact, if U ( ξ ) = 1 M n (cid:80) j ∈ Q nM U j (cid:18) (cid:80) k ∈ Q nM e πi (cid:104) k,j/M − ξ (cid:105) (cid:19) , then U ( (cid:96)/M ) = U (cid:96) ( (cid:96) ∈ Q nM ). It was shown in Section 2.2.4 that if U ( ξ ) jε = m ε ( ξ + v j /
2) (0 ≤ j, ε ≤ n − (1) = 1 and ∂ α m ε ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 = 0 (1 ≤ ε ≤ n − , | α | ≤ d ), then ∂ α m ( ξ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ = v j / = 0 (1 ≤ j ≤ n − , | α | ≤ d ). We conclude that if m (0) = 1 then (cid:88) k ∈ Q nM c αk ( U k ) ε, = ∂ α U ( ξ ) ε, (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 = 0 (1 ≤ ε ≤ n − , | α | ≤ d ) ⇒ ∂ α U ( ξ ) ,j = (cid:88) k ∈ Q nM c αk ( U k ) ,j = 0 (1 ≤ j ≤ n − , | α | ≤ d ) ⇒ (cid:88) k ∈ Q nM c αk U k ∈ C ⊗ C (2 n − × (2 n − ⇒ U ∈ C . We let C (cid:48) = (cid:26) U ∈ ( C n × n ) Q nM σ : (cid:88) k ∈ Q nM c αk U k = (cid:18) a α T c α B α (cid:19) for some a α ∈ C , c α ∈ C n − , B α ∈ C (2 n − × (2 n − (cid:27) . Then we have shown that (cid:18) ∩ n − (cid:96) =0 C ( (cid:96) )1 (cid:19) ∩ C = (cid:18) ∩ n − (cid:96) =0 C ( (cid:96) )1 (cid:19) ∩ C (cid:48) and for this reason, the constraint C may be replaced by C (cid:48) in our algorithms.We now consider the projection onto the subspace described by the regularity constraint C (cid:48) .For k ∈ Q nM , define w k ∈ C n by ( w k ) (cid:96) = ( − (cid:104) k,v (cid:96) (cid:105) , i.e., (cid:8) w k = (1 , ( − (cid:104) k,v (cid:105) , . . . , ( − (cid:104) k,v n − (cid:105) (cid:1) T . We note that because of the definition of the group operation ⊕ on V n defined in (23) we havethat for all k ∈ Z n and integers 0 ≤ j, (cid:96) ≤ n −
1, ( − (cid:104) k,v j ⊕ v (cid:96) (cid:105) = ( − (cid:104) k,v j (cid:105) ( − (cid:104) k,v (cid:96) (cid:105) . Further,from the definition (24) of the permutation matrices σ j , with w k ( k ∈ Q nM ) as above we have( σ j w k ) (cid:96) = n − (cid:88) m =0 ( σ j ) (cid:96)m ( w k ) m = (cid:88) { m : v j ⊕ v (cid:96) = v m } ( − (cid:104) k,v m (cid:105) = ( − (cid:104) k,v j ⊕ v (cid:96) (cid:105) = ( − (cid:104) k,v j (cid:105) ( − (cid:104) k,v (cid:96) (cid:105) = ( − (cid:104) k,v j (cid:105) ( w k ) (cid:96) , so that w k is an eigenvector of σ j with eigenvalue ( − (cid:104) k,v j (cid:105) .We work within the Hilbert space H n = F − M ( C n × n ) Q nM σ = { A ∈ ( C n × n ) Q nM : σ (cid:96) A k = ( − (cid:104) k,v (cid:96) (cid:105) A k for all k ∈ Q nM , 0 ≤ (cid:96) ≤ n − } = (cid:110) A ∈ ( C n × n ) Q nM : A k = w k (cid:0) a k b Tk (cid:1) for some a k ∈ C , b k ∈ C n − and all k ∈ Q nM (cid:111) o that a typical element of H n is a matrix ensemble A ∈ ( C n × n ) Q nM with k -th entry of theform A k = a k b Tk ( − (cid:104) k,v (cid:105) a k ( − (cid:104) k,v (cid:105) b Tk ... ...( − (cid:104) k,v n − (cid:105) a k ( − (cid:104) k,v n − (cid:105) b Tk for some a k ∈ C , b k ∈ C n − . Constraint C (cid:48) is equivalent to the condition (cid:80) k ∈ Q nM k α b Tk = 0( | α | ≤ d ). Let X d be the collection of matrix ensembles B = ( B α ) | α |≤ d ⊂ C n × n of the form B α = (cid:18) γ Tα (cid:19) where = (0 , , . . . , ∈ C n − , ∈ C (2 n − × (2 n − is the zero matrix and γ α ∈ C n − for each | α | ≤ d . We now define an operator R : H n → X d given by( R A ) α = (cid:18) (cid:80) k ∈ Q nM k α b Tk (cid:19) . (47)The projection we require is that onto the kernel of R .We decompose each C ∈ C n × n as C = (cid:18) a b T c D (cid:19) with a ∈ C , b , c ∈ C n − and D ∈ C (2 n − × (2 n − and write a = C , b = C , c = C and D = C . Let A ∈ H n and B ∈ X d .We then have (cid:104)R A , B (cid:105) = (cid:88) | α |≤ d (cid:28)(cid:18) (cid:80) k ∈ Q nM k α b Tk (cid:19) , (cid:18) γ Tα (cid:19)(cid:29) = (cid:88) | α |≤ d (cid:88) k ∈ Q nM k α (cid:10) b Tk , γ Tα (cid:11) = (cid:88) k ∈ Q nM (cid:42) b Tk , (cid:88) | α |≤ d k α γ Tα (cid:43) = (cid:88) k ∈ Q nM (cid:10)(cid:0) a k b Tk (cid:1) , (cid:0) (cid:80) | α |≤ d k α γ Tα (cid:1)(cid:11) = 2 − n (cid:88) k ∈ Q nM (cid:10) w k (cid:0) a k b Tk (cid:1) , w k (cid:0) (cid:80) | α |≤ d k α γ Tα (cid:1)(cid:11) = (cid:104) A , R ∗ B (cid:105) from which we conclude that( R ∗ B ) k = 2 − n w k (cid:0) (cid:80) | α |≤ d k α γ Tα (cid:1) = 2 − n (cid:80) | α |≤ d k α γ Tα − (cid:104) k,v (cid:105) (cid:80) | α |≤ d k α γ Tα − (cid:104) k,v (cid:105) (cid:80) | α |≤ d k α γ Tα ... ...0 ( − (cid:104) k,v n − (cid:105) (cid:80) | α |≤ d k α γ Tα . (48)Let C n,d = |{ α ∈ ( N ∪{ } ) n : | α | ≤ d }| . For example, C ,d = ( d +1)( d +1) /
2. If B = ( B α ) | α |≤ d ∈X d , i.e, B α = (cid:18) γ Tα (cid:19) , then( RR ∗ B ) β = 2 − n (cid:18) (cid:80) | α |≤ d G βα γ Tα (cid:19) here G ∈ C C n,d × C n,d has ( β, α )-th entry G βα = (cid:80) k ∈ Q nM k α + β . We wish to show that RR ∗ isinvertible. Consider functions r α : Q nM → C given by r α ( k ) = k α . We claim that { r α } | α |≤ d isa linearly independent set. To see this, suppose there are constants { a α } | α |≤ d ⊂ C such that (cid:80) | α |≤ d a α r α = 0, i.e., (cid:80) | α |≤ d a α k α = 0 for all k ∈ Q nM . Let p ( x ) = (cid:80) | α |≤ d a α x α ( x ∈ R n ). Then p is a (multivariate) polynomial of degree less than or equal to d and p ( k ) = 0 for all k ∈ Q nM .Hence p ≡
0, i.e., a α = 0 for all α . We conclude that { r α } | α |≤ d is a linearly independent set.Suppose now that a = ( a α ) | α |≤ d is such that G a = . Then0 = ( G a ) β = (cid:88) | α |≤ d G βα a α = (cid:88) k ∈ Q nM k β (cid:88) | α |≤ d a α k α = (cid:28) r β , (cid:88) | α |≤ d a α r α (cid:29) . (49)But (cid:80) | α |≤ d a α r α ∈ S d = sp { r β } | β |≤ d and { r β } | β |≤ d is a basis for S d , so by (49) we conclude that (cid:80) | α |≤ d a α r α ≡
0, or equivalently, p ( k ) = 0 for all k ∈ Q nM where p ( x ) = (cid:80) | α |≤ d a α x α . Hence p ≡ a α = 0 for all α , i.e., G is invertible. We then have(( RR ∗ ) − B ) β = 2 n (cid:18) (cid:80) | α |≤ d G − βα γ Tα (cid:19) . (50)Combining (47), (48) and (50) gives( R ∗ ( RR ∗ ) − R A ) k = (cid:88) | α |≤ d k α (cid:88) | β |≤ d G − αβ (cid:88) (cid:96) ∈ Q nM (cid:96) β (cid:18) b T(cid:96) (cid:19) so that the projection Q of an ensemble A ∈ F − M ( C n × n ) Q nM onto F − M C (cid:48) is given by( Q A ) k = A k − ( R ∗ ( RR ∗ ) − R A ) k = w k (cid:16) a k b Tk − (cid:80) | α |≤ d (cid:80) | β |≤ d (cid:80) (cid:96) ∈ Q nM G − αβ k α (cid:96) β b Tk (cid:17) . Finally, the required projection P C (cid:48) of ( C n × n ) Q nM onto C (cid:48) is given by( P C (cid:48) U ) j = ( F M Q F − M U ) j = U j − (cid:88) k ∈ Q nM ( Q F − M U ) k e πi (cid:104) k,j (cid:105) /M = U j − (cid:88) k ∈ Q nM ( R ∗ ( RR ∗ ) − R ( F − M U )) k e πi (cid:104) k,j (cid:105) /M = U j − M n (cid:88) k ∈ Q nM e πi (cid:104) j,k (cid:105) /M (cid:88) | α |≤ d k α (cid:88) | β |≤ d G − αβ (cid:18) (cid:80) (cid:96) ∈ Q nM (cid:96) β ( F − M U ) T(cid:96), (cid:19) . (51)Let C ∈ C R nd × Q nM have ( β, m )-th entry c βm = (cid:80) (cid:96) ∈ Q nM (cid:96) β e − πi (cid:104) m,(cid:96) (cid:105) /M . Then (51) may be writtenas ( P C (cid:48) U ) j = U j − M n (cid:88) m ∈ Q nM ( C ∗ G − C ) jm (cid:18) U m ) (cid:19) . In this section, we report representative computational results for the DR algorithm (as describedin Algorithm 1) applied to the formulations described in Problem 4.1. The main goal of reporting The accompanying source code is available at https://gitlab.com/matthewktam/drwavelets. hese results is to provide an insight into the typical number of iterations and the success rateof the method for the wavelet reconstruction problem. All experiments implemented in Python3.7 and a machine having an Intel Xeon E5-4650 @ 2.70GHz running Red Hat Enterprise Linux3.10.For each value of ( M, d ) examined, ten replications of the DR algorithm were run, eachstarting from a different randomly generated initialisation x ∈ D , where D denotes the diagonalsubspace from Section 4.2. More precisely, the real and complex entries, respectively, of a matrixensemble U ∈ ( C n × n ) Q nM were generated entry-wise by sampling from the uniform distributionon the interval ( − , +1). The ensemble U was then projected onto C (0)1 , and its projection wasthen used to form the tuple of ensembles x .The algorithm was terminated if either: (i) the stopping criterion (cid:107) x k − x k +1 (cid:107) < (cid:15) was satisfied with (cid:15) = 10 − , or (ii) more than 10 iterations had been performed. In the case thatthe algorithm terminated successfully (i.e., the stopping criterion was satisfied), orthogonalityof the resulting trigonometric polynomial was checked numerically using Bownik’s condition asdescribed in Section 2.2.4. For the 2D problem, non-separability was also checked using theprocedure outlined in Section 2.2.5.Tables 1 and 2 report a summary of the results for the 1D and 2D problems, respectively. Inaddition to the number of instances solved (out of ten), the mean number of iterations and timein seconds for solved instances are shown. The maxima across solved instances are also shown inparentheses. The mean and (in parentheses) maximum separability measure of solved examplesis shown in the final column of table 2.Exemplar results are provided in Figures 1–4. The two-dimensional scaling function andwavelets of Figure 3 and associated filters pass Bownik’s test (Theorem 3.1) for orthogonalityand the separability measure (see Section 2.2.5) of the filter coefficient matrix G is 0 . . H = ( h jk ) j,k =0 satisfying the conditions (cid:88) j,k =0 h jk = 1; (cid:88) j,k =0 | h jk | = 12 . Of course, these filters do not satisfy the extra regularity, consistency, or unitarity conditionssatisfied by the filter given in Figure 3(c). Nevertheless, real-valued scaling functions have beengenerated by the algorithm described in this paper with M = 6, d = 2 and relatively highnon-separability. An example is given in Figure 4. The separability measure of this example isapproximately 0 . Acknowledgement
The authors are grateful for the input of Neil Dizon who helped in the generation of the figuresand provided the highly non-separable example of Section 5.JAH was supported by the Australian Research Council through DP160101537. MKT wassupported by the Australian Research Council through DE200100063. Thanks Roy. ThanksHG. .0 0.2 0.4 0.6 0.8 1.01.000.750.500.250.000.250.500.751.00 ( m )( m )| m | 0.0 0.2 0.4 0.6 0.8 1.01.000.750.500.250.000.250.500.751.00 ( m )( m )| m | (a) The polynomials m ( ξ ) and m ( ξ ), respectively. t t )( t ) (b) The scaling function ϕ and wavelet ψ . [ 0.02490875,-0.0604161 ,-0.09546721,0.3251825 ,0.57055846,0.2352336 ]. (c) The corresponding coefficients. Figure 1: An exemplar 1D result obtained from the DR algorithm for (
M, d ) = (6 , .0 0.2 0.4 0.6 0.8 1.01.000.750.500.250.000.250.500.751.00 ( m )( m )| m | 0.0 0.2 0.4 0.6 0.8 1.01.000.750.500.250.000.250.500.751.00 ( m )( m )| m | (a) The polynomials m ( ξ ) and m ( ξ ), respectively. t t )( t ) 0 1 2 3 4 5 t t )( t ) (b) The scaling function ϕ and wavelet ψ . [-0.046875+0.060515i 0.078125+0.060515i 0.46875 -0.1210307i0.468750-0.121030i 0.078125+0.060515i -0.046875+0.0605153i] (c) The corresponding coefficients. Figure 2: An exemplar 1D result obtained from the DR algorithm for (
M, d ) = (6 , . . . . . . . . . m . . . . . . . . . m . . . . . . . . . m . . . . . . . . . m . . . . . . . ( a ) T h e p o l y n o m i a l s m ( ξ ) , m ( ξ ) , m ( ξ ) a nd m ( ξ ) r e s p ec t i v e l y . T h e r e a l c o m p o n e n t s o f t h e p o l y n o m i a l s a r e s h o w n o n t h e v e r t i c a l a x e s . T h e i m ag i n a r y c o m p o n e n t s a r e r e p r e s e n t e db y c o l o u r . . . . . . . . . . . . ( b ) T h e s c a li n g f un c t i o n ϕ a nd t h e w a v e l e t s ψ ε . T h e r e a l c o m p o - n e n t s o f t h e f un c t i o n s a r e s h o w n o n t h e v e r t i c a l a x e s . T h e i m ag i - n a r y c o m p o n e n t s a r e r e p r e s e n t e db y c o l o u r . [[-1.315e-02+1.780e-02i-2.573e-02+3.276e-02i-1.046e-02+1.251e-02i1.927e-03-2.124e-03i1.794e-04-5.040e-05i3.696e-04-3.775e-04i] [2.368e-02+1.746e-02i4.149e-02+3.268e-02i1.548e-02+1.301e-02i-1.962e-03-1.942e-03i-9.656e-05-2.155e-04i-4.618e-04-4.822e-04i] [1.073e-01-3.551e-02i2.693e-01-6.537e-02i1.573e-01-2.520e-02i-4.749e-02+3.955e-03i-3.024e-02+1.919e-04i1.252e-02+9.006e-04i] [1.070e-01-3.500e-02i2.696e-01-6.551e-02i1.575e-01-2.585e-02i-4.774e-02+4.176e-03i-3.019e-02+3.426e-04i1.249e-02+8.180e-04i] [2.347e-02+1.771e-02i4.167e-02+3.261e-02i1.574e-02+1.269e-02i-2.168e-03-1.831e-03i-1.512e-04-1.415e-04i-4.363e-04-5.232e-04i] [-1.310e-02+1.755e-02i-2.583e-02+3.283e-02i-1.042e-02+1.284e-02i1.965e-03-2.233e-03i7.934e-05-1.271e-04i4.278e-04-3.358e-04i]] ( c ) T h ec o rr e s p o nd i n g s c a li n g f un c t i o n c o e ffi c i e n t s . F i g u r e : A n e x e m p l a r D r e s u l t o b t a i n e d f r o m t h e D R a l go r i t h m f o r ( M , d ) = ( , ) . . . . . . . . . . . m . . . . . . . . . m . . . . . . . . . m . . . . . . . . . m . . . ( a ) T h e p o l y n o m i a l s m ( ξ ) , m ( ξ ) , m ( ξ ) a nd m ( ξ ) r e s p ec t i v e l y . T h e r e a l c o m p o n e n t s o f t h e p o l y n o m i a l s a r e s h o w n o n t h e v e r t i c a l a x e s . T h e i m ag i n a r y c o m p o n e n t s a r e r e p r e s e n t e db y c o l o u r . . . . . . . . . . . . . . . ( b ) T h e ( r e a l - v a l u e d ) s c a li n g f un c t i o n ϕ . [[0.0368,0.0406,-0.0305,-0.0362,0.0061,0.008], [-0.0341,-0.0591,-0.0079,0.0292,0.0118,-0.0003], [-0.0308,-0.0405,-0.0225,-0.0185,0.0056,0.0113], [0.0304,0.0661,0.082,0.0877,0.0502,0.0088], [0.0065,-0.0303,0.0052,0.2172,0.2736,0.0984], [0.0161,-0.0371,-0.1218,0.0457,0.2233,0.1091]] ( c ) T h ec o rr e s p o nd i n g s c a li n g f un c t i o n c o e ffi c i e n t s . F i g u r e : A n e x e m p l a r D r e s u l t o b t a i n e d f r o m t h e D R a l go r i t h m f o r ( M , d ) = ( , ) . (cid:15) = 10 − .( M, d ) Solved Iterations Time (s)(4 ,
1) 10 122.2 (162) 0.1 (0.2)(6 ,
2) 9 3 852.0 (9 361) 5.1 (12.5)(8 ,
3) 10 40 672.5 (112 460) 67.6 (186.7)(10 ,
4) 8 154 372.8 (607 495) 325.8 (1 280.5)(12 ,
5) 9 166 251.0 (369 136) 422.0 (932.8)(14 ,
6) 6 302 014.3 (690 650) 917.1 (2 093.7)Table 2: Mean (worst case) results from 10 replications for the 2D problem with (cid:15) = 10 − .( M, d ) Solved Iterations Time (s) S ( ϕ )(4 ,
1) 10 4 469.2 (28 387) 118.9 (708.2) 0.209 (0.250)(6 ,
2) 6 180 864.3 (747 870) 22 322.2 (92 288.1) 0.104 (0.207)
References [1] Arag´on Artacho F. J., Borwein i. M. & Tam M. K. (2013) Recent results on Douglas–Rachford methods for combinatorial optimization problems,
Journal of Optimization Theoryand Applications , 163(1):1–30.[2] Arag´on Artacho F. J., Borwein i. M. & Tam M. K. (2014) Douglas–Rachford feasibilitymethods for matrix completion problems,
The ANZIAM Journal , 55(4):299–326.[3] Arag´on Artacho F. J. & Campoy R. (2018) A new projection method for finding the closestpoint in the intersection of convex sets,
Computational Optimization and Applications ,69(1):99–132.[4] Arag´on Artacho F. J., Campoy R. & Elser V. (2020) An enhanced formulation for solvinggraph coloring problems with the Douglas–Rachford algorithm,
Journal of Global Optimiza-tion .[5] Arag´on Artacho F. J., Campoy R., Kostsireas I. & Tam M. K. (2018) A feasibility ap-proach for constructing combinatorial designs of circulant type,
Journal of CombinatorialOptimization , 35(4):1061–1085.[6] Ayache A. (1999) Construction of non separable dyadic compactly supported orthonormalwavelet bases for L ( R ) of arbitrarily high regularity, Revista Matem´atica Iberoamericana ,15(1):37–58[7] Bauschke H. H. & Combettes P. L. (2011)
Convex analysis and monotone operator theoryin Hilbert spaces , New York: Springer.[8] Bauschke H. H. Koch V. R. & Phan H. M. (2016) Stadium Norm and Douglas–RachfordSplitting: A New Approach to Road Design Optimization,
Operations Research , 64(1):201–218.[9] Bauschke H. H., Combettes P. L. & Luke D. R. (2004)Finding best approximation pairs relative to two closed convex sets in Hilbert spaces,
Journal of Approximation Theory , 127(2):178–192.
10] Belogay B. & Wang Y. (1999) Arbitrarily Smooth Orthogonal Nonseparable Wavelets in R , SIAM Journal on Mathematical Analysis , 30(3):678–697.[11] Borwein J. M. & Tam M. K. (2014) A cyclic Douglas–Rachford iteration scheme,
Journalof Optimization Theory and Applications , 160:1–29.[12] Bregman L. M. (1965) The method of successive projection for finding a common point ofconvex sets,
Doklady Akademii Nauk , 162(3):688–692.[13] Bownik M. (1997) Tight frames of multidimensional wavelets,
Journal of Fourier Analysisand Applications . 3(5):525-542.[14] Calder´on A.P. (1964) Intermediate spaces and interpolation, the complex method,
StudiaMathematica , 24:113–190.[15] Cohen A. (1990) Ondelettes, analys´ees multir´esolutions et filtres miroir en quadrature,
Ann.Inst. H. Poincar´e, Anal non lin´eaire
RevistaMatematica Iberoamericana . 9(1):51–137.[17] Dao M. N. & Tam M. K. (2019) Union averaged operators with applications to proximalalgorithms for min-convex functions,
Journal of Optimization Theory and Applications
Comm. PureAppl. Math.
Ten Lectures on Wavelets,
SIAM.[20] Dizon N.D., Hogan J.A. & Lakey J.D. (2019) Optimization in the construction of nearlycardinal and nearly symmetric wavelets, 2019 International Conference on Sampling Theoryand Applications (SampTA), Bordeaux, France[21] Duffin R.J. & Schaeffer A.C. (1952) A class of nonharmonic Fourier series,
Transactions ofthe American Mathematical Society
72: 341–366.[22] Gravel S. & Elser V. (2008) Divide and concur: A general approach to constraint satisfac-tion,
Physical Review E , 78(3):036706.[23] Franklin D.J. (2018)
Projection Algorithms for Non-Separable Wavelets and Clifford FourierAnalysis,
PhD Thesis, University of Newcastle, Australia.[24] Franklin D.J., Hogan J.A. & Tam, M. (2019) Higher dimensional wavelets and the Douglas–Rachford algorithm, 2019 International Conference on Sampling Theory and Applications(SampTA), Bordeaux, France.[25] Gilbert J.E., Han Y.S., Hogan J.A., Lakey J.D., Weiland D. & Weiss G. (2002).
SmoothMolecular Decompositions of Functions and Singular Integral Operators,
Memoirs of theAmerican Mathematical Society 156 (1). Newport, Rhode Island: American MathematicalSociety.[26] Grossman J.& Morlet J. (1985) Decompositions of Hardy functions into square integrablewavelets of constant shape,
SIAM journal on Mathematical Analysis
IEEE Transactions on Image Processing , 9(5):949–953.
28] Hesse R., Luke D.R. & Neumann P. (2014). Alternating projections and Douglas–Rachfordfor sparse affine feasibility,
IEEE Transactions on Signal Processing , 62(18):4868–4881.[29] Higham N.J. (2008)
Functions of matrices: theory and computation,
SIAM.[30] Hogan J.A. & Lakey J.D. (2005)
Time-Frequency and Time-Scale Methods : Adaptive De-compositions, Uncertainty Principles, and Sampling,
John Benedetto (Ed.), Basel, Switzer-land: Birkh¨auser.[31] Karoui A. (2003) A note on the construction of nonseparable wavelet bases and multiwaveletmatrix filters of L ( R n ), where n ≥ Electronic Research Announcements of the AmericanMathematical Society. L ( R ), Applied Mathematics Letters. R n , IEEE Transactions on Information Theory , 38(2):533-555.[34] Lai M. & Roach D.W. (1999) Nonseparable symmetric wavelets with short support, in
Wavelet Applications in Signal and Image Processing VII . International Society for Opticsand Photonics, 3813:132–147.[35] Lai M.-i. (2002) Methods for constructing nonseparable compactly supported orthonormalwavelets, in
Wavelet Analysis: Twenty Years’ Developments.
SIAM Journal on Numerical Analysis , 16(6):964–979.[37] Mallat S. (1989) Multiresolution approximation and wavelets,
Transactions of the AmericanMathematical Society , 315:69–88.[38] Meyer Y. (1986) Ondelettes, fonctions splines et analyses gradu´ees, Lectures given at theUniversity of Torino.[39] Meyer Y. (1989)
Wavelets and Operators,
Cambridge University Press.[40] Phan H.M. (2016). Linear convergence of the Douglas–Rachford method for two closed sets,
Optimization , 65(2):369–385.[41] San Antolin A. & Zalik R.A. (2013) A family of nonseparable scaling functions and com-pactly supported tight framelets,
Journal of Mathematical Analysis and Applications.
Introduction to Fourier Analysis on Euclidean Spaces,