A General Approach for Orthogonal 4-Tap Integer Multiwavelet Transforms

Abstract

An algorithm for orthogonal 4-tap integer multiwavelet transforms is proposed. We compute the singular value decomposition (SVD) of block recursive matrices of transform matrix, and then transform matrix can be rewritten in a product of two block diagonal matrices and a permutation matrix. Furthermore, we factorize the block matrix of block diagonal matrices into triangular elementary reversible matrices (TERMs), which map integers to integers by rounding arithmetic. The cost of factorizing block matrix into TERMs does not increase with the increase of the dimension of transform matrix, and the proposed algorithm is in-place calculation and without allocating auxiliary memory. Examples of integer multiwavelet transform using DGHM and CL are given, which verify that the proposed algorithm is an executable algorithm and outperforms the existing algorithm for orthogonal 4-tap integer multiwavelet transform.

1. Introduction

In many applications of image processing, the given data are integer valued. To compress digital image losslessly by means of transformation, the transform must map integers to integers and be perfectly invertible. Calderbank et al. [1] presented two wavelet transforms approaches that map integers to integers, and the proposed lifting scheme has been a popular method for mapping integers to integers. Using lifting scheme, Deever and Hemami [2] presented a projection-based technique for decreasing the first-order entropy of transform coefficients. However, Multiwavelets have several advantages in comparison to scalar wavelets. Such features as short support, orthogonality, symmetry, and vanishing moments are known to be important in signal processing. A scalar wavelet cannot possess all these properties at the time.

Unlike scalar wavelet, little literature is available on integer multiwavelet transform. Cheung et al. [3] presented an integer multiwavelet transform and its associated integer prefilter based on box-and–slope multiscaling system. Van Fleet [4] developed a factorization of 4-tap multiwavelets with multiplicity 2. In particular, the author in [4] derived a factorization of the DGHM multiwavelet and applied DGHM multiwavelet processing to data compression. In order to implement the DGHM integer transform, the author computed scaling value 𝛼 under certain conditions [4].The scaling value 𝛼 is not convenient to compute and larger than 1, which turns out to be bad for lossless compression [1].

This paper is to present an algorithm of 4-tap orthogonal multiwavelets that map integers to integers and be perfectly invertible using SVD and TERMs. Unlike [4], the proposed algorithm does not compute scaling value 𝛼 and calculates in-place and without allocating auxiliary memory.

2. Integer Transform and Multiwavelet

2.1. Integer Transform

Hao and Shi [5] proved that there exists an implementation of integer mapping when a linear transform is invertible and in finite-dimensional space, as follows.

Theorem 2.1. Matrix 𝐴 has a TERM factorization of 𝐴=𝑃𝑉₁𝑉₂⋯𝑉_𝑀𝐷_𝑅 if and only if |det(𝐴)|=1, where 𝑀 is finite, 𝑉_𝑘(𝑘=1,2…𝑀) are unit TERMs, 𝑃 is a permutation matrix, and 𝐷_𝑅 is a rotator for only one complex number. If all the diagonal elements of a TERM are equal to 1, the TERM will be a unit triangular matrix. If 𝐴 is a real matrix, then 𝐷_𝑅 is an identity matrix.

Unless otherwise stated, we will concentrate on the nonsingular real matrix 𝐴, thus reversible integer implementation for general linear transforms 𝑌=𝐴𝑋 can be derived as follows [5], where 𝑌=[𝑦₁,𝑦₂,…,𝑦_𝑛]^T,  𝑋=[𝑥₁,𝑥₂,…,𝑥_𝑛]^T.

Let 𝐿=(𝑙_𝑖,𝑗)_𝑛×𝑛,𝑈=(𝑢_𝑖,𝑗)_𝑛×𝑛, and 𝑃_𝑛×𝑛 be a lower TERM, an upper TERM, and a permutation matrix, respectively. If 𝐴=𝐿 is a lower TERM, the computational ordering of linear transform 𝑌=𝐴𝑋 can be arranged to be top-down:𝑦_𝑖=𝑙_𝑖,𝑖𝑥_𝑖𝑙,𝑖=1,_𝑖,𝑖𝑥_𝑖+[^𝑖−1_𝑘=1𝑙_𝑖,𝑘𝑥_𝑘],2≤𝑖≤𝑛. (2.1) Its inverse ordering is reversed:𝑥_𝑖=𝑦_𝑖𝑙_𝑖,𝑖(1,𝑖=1,𝑙_𝑖,𝑖)(𝑦_𝑖−[^𝑖−1_𝑘=1𝑙_𝑖,𝑘𝑥_𝑘]),2≤𝑖≤𝑛. (2.2)

Likewise, if 𝐴=𝑈 is an upper TERM, the computational ordering of linear transform 𝑌=𝐴𝑋 can be arranged to be top-down:𝑦_𝑖=𝑢_𝑖,𝑖𝑥_𝑖𝑢,𝑖=𝑛,_𝑖,𝑖𝑥_𝑖+[^𝑛_𝑘=𝑖+1𝑢_𝑖,𝑘𝑥_𝑘],1≤𝑖≤𝑛−1. (2.3) Its inverse ordering is reversed:𝑥_𝑖=𝑦_𝑖𝑢_𝑖,𝑖(1,𝑖=𝑛,𝑢_𝑖,𝑖)(𝑦_𝑖−[^𝑛_𝑘=𝑖+1𝑢_𝑖,𝑘𝑥_𝑘]),1≤𝑖≤𝑛−1, (2.4) where [] denotes rounding arithmetic. Since permutation matrix 𝑃_𝑛×𝑛 is an orthogonal matrix, transform and the inverse transform can be denoted by 𝑌=𝑃_𝑛×𝑛𝑋 and 𝑋=𝑃^T_𝑛×𝑛𝑌, respectively.

2.2. Properties of Transform Matrix

While many researchers have investigated multiwavelet of multiplicity 𝑟 [6–9], little work has been published on integer multiwavelet transforms. We will discuss integer transform on 4-tap multiwavelet of multiplicity 𝑟. Let Φ(𝑡)=[𝜙₁,…𝜙_𝑟]^T,  Ψ(𝑡)=[𝜓₁,…𝜓_𝑟]^T∈𝐿²(𝑅)^𝑟 be a multiscaling function and a multiwavelet function, respectively, and then both Φ(𝑡) and Ψ(𝑡) satisfy the following two-scale relation: Φ(𝑡)=³_𝑘=0𝑃_𝑘Φ(2𝑡−𝑘),Ψ(𝑡)=³_𝑘=0𝑄_𝑘Φ(2𝑡−𝑘), (2.5) where 𝑃_𝑘 and 𝑄_𝑘 are 𝑟×𝑟 matrices, respectively. For more details about multiwavelets, see, for example, [7–9]. As is the case in [4], the matrix representation 𝐻_𝑛×𝑛 of multiwavelet transformation is𝐻_𝑛×𝑛=𝐻₀𝐻₁𝐻₀𝐻₁𝐻⋱⋱₀𝐻₁𝐻₁𝐻₀, (2.6) where 𝐻_𝑖=[^{𝑃2𝑖𝑃2𝑖+1}_{𝑄2𝑖𝑄2𝑖+1}],  𝑖=0,1,𝑃_𝑗,𝑄_𝑗,𝑗=0,1,2,3, are 𝑟×𝑟 matrices, respectively.

In order to find integer multiwavelet transform schemes, inspired by [1], we will prove some remarkable properties of orthogonal matrix 𝐻 and then present the approach that maps integers to integers.

Property 1. Suppose that rank(𝐻₀)=𝑚, rank(𝐻₁)=𝑠, thus 𝑚+𝑠=2𝑟.

Proof. From the condition of orthogonality, 𝐻^T𝐻=𝐻𝐻^T=𝐼_𝑛×𝑛. It follows that 𝐻₀𝐻^T₀+𝐻₁𝐻^T₁=𝐼_2𝑟×2𝑟,𝐻₁𝐻^T₀=𝐻₀𝐻^T₁=0,2𝑟=rank(𝐼)=rank(𝐻₀𝐻^T₀+𝐻₁𝐻^T₁)≤rank(𝐻₀)+rank(𝐻₁),0=rank(𝐻₀𝐻^T₁)≥rank(𝐻₀)+rank(𝐻₁)−2𝑟. (2.7) Thus, rank(𝐻₁)+rank(𝐻₀)=2𝑟. The proof is completed.

Let 𝑈₀Σ₀𝑉^T₀ and 𝑈₁Σ₁𝑉^T₁ be the singular value decomposition of 𝐻₀,𝐻₁, respectively. By Property 1, since rank(𝐻₁)=rank(Σ₁), rank(𝐻₀)=rank(Σ₀), we see that rank(Σ₀)+rank(Σ₁)=2𝑟.

Property 2. If the 𝜆_𝑖 are the singular values of 𝐻₀ and the vectors 𝑢^𝑖₀ and 𝑣^𝑖₀ are the 𝑖th left singular vector and 𝑖th right singular vector, respectively, and the 𝜇_𝑖 are the singular values of 𝐻₁ and the vectors 𝑢^𝑗₁ and 𝑣^𝑗₁ are the 𝑗th left singular vector and 𝑖th right singular vector, respectively, then 𝑢^𝑖₀ and 𝑣^𝑖₀ are orthogonal to the vectors 𝑢^𝑗₁ and 𝑣^𝑗₁, respectively, 𝑖=1,2,…𝑚,  𝑚=rank(𝐻₀),  𝑗=1,2…𝑠,𝑠=rank(𝐻₁).

Proof. From (2.7), we have 𝐻₁𝐻^T₀=𝑈₁Σ₁𝑉^T₁𝑉₀Σ₀𝑈^T₀=0. (2.8) Since 𝑈₀,𝑈₁ are two orthogonal matrices, we have Σ₁𝑉^T₁𝑉₀Σ₀=0,Σ₁𝑉^𝑇₁𝑉₀Σ₀=Σ₁[𝑣¹₁,…,𝑣₁^2𝑟]^𝑇[𝑣¹₀,…,𝑣₀^2𝑟]Σ₀=𝜆₁⋱𝜆_𝑚0(𝑣¹₁)^T(𝑣²₁)^T⋮(𝑣₁^2𝑟)^T𝑣¹₀𝑣²₀,…,𝑣₀^2𝑟𝜇₁⋱𝜇_𝑠0=𝜆₁(𝑣¹₁)^T⋮𝜆_𝑚(𝑣^𝑚₁)^T0𝜇₁𝑣¹₀,…,𝜇_𝑠𝑣^𝑠₀0=𝜆₁𝜇₁⟨𝑣¹₁,𝑣¹₀⟩⋯𝜆₁𝜇_𝑠⟨𝑣¹₁,𝑣^𝑠₀𝜆⟩0⋮⋮⋮_𝑚𝜇₁⟨𝑣^𝑚₁,𝑣¹₀⟩⋯𝜆_𝑚𝜇_𝑠⟨𝑣^𝑚₁,𝑣^𝑠₀⟩00⋯00=0. (2.9)(1)(2)(3)(4) From (2.9), we have ⟨𝑣^𝑖₀,𝑣^𝑗₁⟩=0. Likewise, ⟨𝑢^𝑖₀,𝑢^𝑗₁⟩=0. The proof is completed.

By Property 2, we can rewrite the matrices 𝐻₀ and 𝐻₁ in this case as𝐻₀=𝑈₀Σ₀𝑉^T₀=[𝑢¹₀,𝑢²₀,…,𝑢^𝑚₀,𝑢¹₁,𝑢²₁,…,𝑢^𝑠₁]×diag(𝜆₁,…,𝜆_𝑚,0⋯0_𝑠)[𝑣¹₀,𝑣²₀,…,𝑣^𝑚₀,𝑣¹₁,𝑣²₁,…,𝑣^𝑠₁]^T=𝑈Σ₀𝑉^T,𝐻₁=𝑈₁Σ₁𝑉^T₁=[𝑢¹₀,𝑢²₀,…,𝑢^𝑚₀,𝑢¹₁,𝑢²₁,…,𝑢^𝑠₁]diag(0,…,0_𝑚,𝜇₁,…,𝜇_𝑠)×[𝑣¹₀,𝑣²₀,…,𝑣^𝑚₀,𝑣¹₁,𝑣²₁,…,𝑣^𝑠₁]^T=𝑈Σ₁𝑉^T, (2.10)(2.11) where 𝑈=[𝑢¹₀,𝑢²₀,…,𝑢^𝑚₀,𝑢¹₁,𝑢²₁,…,𝑢^𝑠₁],  𝑉=[𝑣¹₀,𝑣²₀,…,𝑣^𝑚₀,𝑣¹₁,𝑣²₁,…,𝑣^𝑠₁]∑₀=diag(𝜆₁,…,𝜆_𝑚,0⋯0_𝑠)  ∑₁=diag(0⋯0,_𝑚𝜇₁,…,𝜇_𝑠). From Property 2, we see that 𝑈 and 𝑉 are two orthogonal matrices.

Property 3. If the 𝜆_𝑖 and 𝜇_𝑗 are the nonzero singular values of 𝐻₀ and 𝐻₁, respectively, then the nonzero singular values of 𝐻₀ and 𝐻₁ are 1.

Proof. From (2.7), we can derive 𝐼=𝐻₀𝐻^T₀+𝐻₁𝐻^T₁=𝑈Σ₀𝑉^T𝑉Σ₀𝑈^T+𝑈Σ₁𝑉^T𝑉Σ₁𝑈^T=𝑈(Σ²₀+Σ²₁)𝑈^T. (2.12) Left multiplying (2.12) with 𝑈⁻¹and right multiplying (2.11) with(𝑈^T)⁻¹, we have 𝐼=Σ²₀+Σ²₁=diag(𝜆²₁,…,𝜆²_𝑚,𝜇²₁,…,𝜇²_𝑠). The proof is completed.

By Properties 1 and 2, we can rewrite the matrix 𝐻_𝑛×𝑛 in this case as𝐻_𝑛×𝑛=𝐻₀𝐻₁𝐻₀𝐻₁⋱𝐻₁𝐻₀=𝑈𝑈⋱𝑈Σ₀Σ₁Σ₀Σ₁Σ₁Σ₀𝑉^T𝑉^T⋱𝑉^T=𝑈𝑃𝑉^T, (2.13) where 𝑈_𝑛×𝑛 and 𝑉^T_𝑛×𝑛 are block diagonal with the same 2𝑟×2𝑟 orthogonal matrices 𝑈,𝑉^T comprising the diagonal, respectively, and 𝑃_𝑛×𝑛 is a permutation matrix. The advantages of the re-representation of transform 𝐻_𝑛×𝑛 are as follows:(1)the TERMs factorization of 𝐻_𝑛×𝑛 directly requires more calculation than 𝑈 and 𝑉, because the TERMs factorization of 𝐻_𝑛×𝑛 directly involves 𝑜(𝑛²) operations [5], and the TERMs factorization of 𝑈(𝑉) involves 𝑜((2𝑟)²) operations; in general, 𝑈(𝑉) is 4 by 4 matrices with the case 𝑟=2;(2)integer multiwavelet transform can map integers to integers in-place and without allocating auxiliary memory (see Section 2.1).

3. Multiwavelet Integer Transform

Since 𝑈,𝑉 are two 2𝑟×2𝑟 orthogonal matrices, we suppose that the input data {𝑥₁,𝑥₂,…,𝑥_𝑛} is denoted by {[𝑥₁^1,𝑗,𝑥₂^1,𝑗,…𝑥^1,𝑗_2𝑟]^T,[𝑥₁^2,𝑗,𝑥₂^2,𝑗,…𝑥^2,𝑗_2𝑟]^T,…,[𝑥₁^𝐿,𝑗,𝑥₂^𝐿,𝑗,…𝑥^𝐿,𝑗_2𝑟]^T}, that is, {𝑋^1,𝑗,𝑋^2,𝑗,…𝑋^𝐿,𝑗}, where 𝐿 is a positive integer, the superscript 𝑗 denotes the 𝑗th stages of the transform of 𝑋^𝑘,𝑗, 𝑘=1,2,…𝐿, and 𝑘 is the 𝑘th sequence of input data. When {𝑋^𝑘,𝑗}^𝐿_𝑘=1 is an initial input data, 𝑗 equals 0. We can state an algorithm for 4-tap orthogonal integer multiwavelet transform as follows.(1)Compute the singular value decomposition of 𝐻₀, 𝐻₁, and find 𝑈,  𝑃=[Σ₀|Σ₁]_2𝑟×4𝑟,  𝑉^T.(2)Find the factorization of 𝑈, 𝑉^T, such as 𝑈=𝑃₀𝑉₁𝑉₂,…,𝑉_𝑀, 𝑉^T=𝑃₁𝑊₁𝑊₂,…,𝑊_𝑀, respectively.(3)Combine the data sequence {𝑋^𝑘,𝑗}^𝐿_𝑘=1with the matrices 𝑊_𝑀,…𝑊₁,𝑃₁,𝑃,𝑉_𝑀,…𝑉₁,𝑃₀ by corresponding equations in Section 2.1, sequentially, according to transform matrix form (upper, lower, and permutation matrix).

Specially speaking, when 𝑗 equals 𝑀+1, 𝑃 is a 2𝑟×4𝑟 matrix, so the ordering of transform can be carried through𝑋^𝑘,𝑗+1=Σ₀𝑋^𝑘,𝑗+Σ₁𝑋^𝑘+1,𝑗Σ,𝑘=1,2…(𝐿−1),₀𝑋^𝑘,𝑗+Σ₁𝑋^1,𝑗,𝑘=𝐿, (3.1) where 𝑋^𝑘,𝑗 is a 2𝑟×1 matrix.

The 𝑗th stage of the transform is a product of 𝑋^𝑘,𝑗 and the 𝑗th matrix of the matrices 𝑊_𝑀,…𝑊₁,𝑃₁,𝑃,𝑉_𝑀,…𝑉₁,𝑃₀ from left to right by rounding arithmetic (see Section 2.1). Since 𝑋^𝑘,𝑗 is integer valued, 𝑋^𝑘,𝑗+1 is also integer valued. After 𝑗=𝑀+𝑀+3 stage of the transform, each stage of this process can map integers to integers, therefore transform results of the initial data sequence {𝑋^𝑘,𝑗}^𝐿_𝑘=1 are integer valued.

The inverse transform must be the backward running of the forward transformation.

4. Experimental Results

In order to demonstrate the availability of the proposed scheme, DGHM [6] and CL [7] multiwavelets are used as examples. The preprocessing is given by balancing the nonbalanced multiwavelet. The way to achieve this is to find the orthogonal matrix 𝑅, such that 𝑅^T𝑃_𝑖𝑅=𝑃_𝑖 and 𝑅^T𝑄_𝑖𝑅=𝑄_𝑖 [10]. Since 𝑅 is an orthogonal matrix, the new matrices 𝐻₀=𝐷^T𝐻₀𝐷 and 𝐻₁=𝐷^T𝐻₁𝐷 inherit all of the properties of 𝐻₀ and 𝐻₁ in Section 2.2, where 𝐷=diag(𝑅,𝑅). We can compute the SVD of matrices 𝐻₀ and 𝐻₁, respectively, thus we will have two matrices 𝑈 and 𝑉^T.

For the DGHM case, the corresponding matrices in (2.5) are 𝑃₀=3√24105−1√20−3220,𝑃₁=900√20−32,𝑃20₂=3√20910√2022,𝑃₃=00−10,𝑄20₀=√33√2062000,𝑄₁=−9√3√2066√0−33,𝑄₂=3√3−√2063√206−√1035,𝑄₃=−√30−√6060.30 (4.1) By orthogonal matrix 𝑅 [11], 1𝑅=√61−√21+√21+√2−1+√2, (4.2) we can find the matrices 𝑈,  𝑉^T, and 𝑃 as follows: 𝑈=,0.01579120.1188553−0.989179−0.084551−0.6383476−0.5786332−0.1218390.49279930.7295663−0.6631722−0.08055460.14644660.24494070.4596299−0.0138210.8535535𝑉^T=,00.3445041−0.8684815−0.35645020−0.07195550.3541483−0.932417−0.2841407−0.8974428−0.3325729−0.05706050.9587826−0.2659623−0.0985599−0.0169102 (4.3)𝑃=[Σ₀|Σ₁], where Σ₀=diag(1,1,1,0), Σ₁=diag(0,0,0,1), and then compute the TERMs factorization of 𝑈,  𝑉^T as follows: 𝑈=𝑃₀𝑉₁𝑉₂,…,𝑉₈=×××00100100100000011000−0.5042912100−0.0209632−0.1149459100.1175034−0.58873660.2706307111000100001100011000−0.2704337100001000−0.146446711−1.3706774000100001−1.17157280001100001100−2.18275211000011000011000100001100001000−1.84548571000011−0.9089951−0.11041440.2007310011.3487059−0.51665130010.01619230001,𝑉^T=𝑃₁𝑊₁𝑊₂,…,𝑊₈=×××00100001010010001000−0.2533664100−0.0209632−0.3528807100.00316850.0737051−0.3098048111000100001100011000−0.0412174100001000−0.06758311−1.0429893000100001−1.07248150001100001000−2.06834941000011000011000100001100001000−1.93602321000011−0.2773958−0.1027969−0.0176372011.43892790.46935910010.37981750001. (4.4)(5)(6)(7)(8)(9)(10)(11)(12)(13) Likewise, for the CL case, the corresponding matrices in (2.5) are𝑃₀=10−3√105√406−2√155√406−3√15√405−31040,𝑃₁=30+3√105√406−2√15−5√406+7√15√4015−310𝑃40₂,=30+3√10−5√406−2√155√406+7√15√4015−31040,𝑃₃=10−3√10−5√406−2√15−5√406−3√15√405−310𝑄40₀,=5√6−2√15−√4010−310−√405−3105√406−3√1540,𝑄₁=−5√6−2√15√4030+310√4015−3105√406+7√15𝑄40₂,=−5√6−2√15−√4030+310−√4015−3105√406+7√1540,𝑄₃=5√6−2√15√4010−310√405−3105√406−3√1540. (4.5)(14)(15)(16) By matrix 𝑅 [11],√𝑅=221−111, (4.6) matrices 𝑈,  𝑉^𝑇, and 𝑃 can be represented as𝑈=,−0.7104132−0.6924131−0.09023430.0879480.087948−0.0902343−0.6924132−0.71041330.09023430.087948−0.71041320.69241320.6924132−0.71041330.0879480.0902343𝑉^T=,0.11326590−0.993564700−0.11326590−0.9935648−0.9935647−0.0000001−0.11326590.00000010.0000001−0.993564700.1132659 (4.7)𝑃=[Σ₀|Σ₁], where Σ₀=diag(1,1,0,0), Σ₁=diag(0,0,1,1), and then compute the TERMs factorization of 𝑈, 𝑉^T as follows: 𝑈=𝑃₀𝑉₁𝑉₂⋯𝑉₈=××10000001001001001000−3.3822942100−0.1270166010−0.42960760.1270166−1.4329687111000100001100011000−1.7104132100001000−1.7104132111.40763150001000011.407631400011000010000.01613321000011000011000100001100001000−0.01587710000110.97466250.1270166−0.123798401−1.01613320.8633704001−0.97466270001,𝑉^T=𝑃₁𝑊₁𝑊₂,…,𝑊₈=××00100001100001001000−2.0064769100−0.1139995010−0.22873740.1139995−2.0064769111000100001100011000−1.9935647100001000−1.9935648111.0064770001000011.006476900011000011000100001100001000−0.0128291100001100.1139995001−1.012996−0.1139994001−0.000000100011000010000.0129958100001. (4.8)(17)(18)(19)(20)(21)(22)(23) In order to test the effectiveness of the proposed scheme, the effectiveness for lossless compression is measured using the weighted entropy given by [1]. For comparison, we will focus on 4-tap multiwavelets, such as CL, DGHM, and two scalar wavelets: the Daubechies four-tap orthogonal wavelet D4 and the Daubechies biorthogonal wavelet 9-7. We use seven test images, which are 512×512 pixels in size, and we decompose each image one level with each transform. The weighted entropies are tabulated in Table 1.

Table 1 

Weighted entropies of transformed images.

These results suggest that our method outperforms the method of [4], because the method of [4] enlarges the dynamic range of the output data by scaling value 𝛼. It is also worth noting that the performance of CL is very close to that of 9-7. Since preprocessing of the input data is a crucial point in nonbalanced multiwavelet applications [8], we used order-one balanced multiwavelets by balancing the nonbalanced multiwavelet [10]. The focus of this paper is to develop an algorithm of integer multiwavelet transform rather than construct several balanced multiwavelets by orthogonal matrix 𝑅, and then a further improvement in the balanced order is needed.

5. Conclusions

This paper presented a novel algorithm for orthogonal 4-tap integer multiwavelet transforms. Unlike the existing algorithm, transform matrix can be rewritten in a product of two block diagonal matrices and a permutation matrix, and then we factorize block matrix of block diagonal matrices into TERMs, which is lower than the TERMs factorization of 𝐻_𝑛×𝑛 directly in computational cost. In addition, the proposed algorithm computes in-place and without allocating auxiliary memory. Better results could perhaps be achieved with higher balanced order and also with other multiwavelet bases.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 60972142) and 973 program (Project no. 2010CB327904).