On the Low-Rank Approximation Arising in the Generalized Karhunen-Loeve Transform

Abstract

We consider the low-rank approximation problem arising in the generalized Karhunen-Loeve transform. A sufficient condition for the existence of a solution is derived, and the analytical expression of the solution is given. A numerical algorithm is proposed to compute the solution. The new algorithm is illustrated by numerical experiments.

1. Introduction

Throughout this paper, we use 𝑅^𝑚×𝑛 to denote the set of 𝑚×𝑛 real matrices. We use 𝐴^𝑇 and 𝐴⁺ to denote the transpose and Moore-Penrose generalized inverse of the matrix 𝐴, respectively. The symbol 𝑂^𝑛×𝑛 stands for the set of all 𝑛×𝑛 orthogonal matrices. The symbols rank(𝐴) and ‖𝐴‖_𝐹 stand for the rank and the Frobenius norm of the matrix 𝐴, respectively. For 𝑎=(𝑎_𝑖)∈𝑅^𝑛, the symbol ‖𝑎‖ stands for the 𝑙₂-norm of the vector 𝑎, that is, ‖𝑎‖₂=(∑^𝑛_𝑖=1𝑎²_𝑖)^1/2. The symbol 𝐴^1/2 stands for the square root of the matrix 𝐴, that is, (𝐴^1/2)²=𝐴. For the random vector 𝑥=(𝑥_𝑖)∈𝑅^𝑛, we use 𝐸{𝑥_𝑖} to stand for the expected value of the 𝑖th entry 𝑥_𝑖, and we use 𝐸{𝑥𝑥^𝑇}=(𝑒_𝑖𝑗)_𝑛×𝑛 to stand for the covariance matrix of the random vector 𝑥, where 𝑒_𝑖𝑗=𝐸[(𝑥_𝑖−𝐸{𝑥_𝑖})(𝑥_𝑗−𝐸{𝑥_𝑗})], 𝑖,𝑗=1,2,…,𝑛.

The generalized Karhunen-Loeve transform is a well-known signal processing technique for data compression and filtering (see [1–4] for more details). A simple description of the generalized Karhunen-Loeve transform is as follows. Given two random vectors 𝑥∈𝑅^𝑛, 𝑠∈𝑅^𝑚 and an integer 𝑑 (1≤𝑑<min{𝑚,𝑛}), the generalized Karhunen-Loeve transform is presented by a matrix 𝑇^∗, which is a solution of the following minimization problem (see [1, 4]): min_{𝑇∈𝑅𝑚×𝑛(𝑇),rank=𝑑}𝐸‖𝑠−𝑇𝑥‖². (1) Here the vector 𝑠 depends on some prior knowledge about the data 𝑥.

Without the rank constraint on 𝑇, the solution of the minimization problem (1) is 𝑇₀=𝑅_𝑠𝑥𝑅⁺_𝑥, (2) where 𝑅_𝑠𝑥=𝐸{𝑠𝑥^𝑇}, 𝑅_𝑥=𝐸{𝑥𝑥^𝑇}. The minimization problem with this case is associated with the well-known concept of Wiener filtering (see [3]).

With the rank constraint on 𝑇, that is, rank(𝑇)=𝑑, we first consider the cost function of the minimization problem (1). By using the fact 𝑅_𝑠𝑥𝑅_𝑥𝑅⁺_𝑥=𝑅_𝑠𝑥 and the four Moore-Penrose equations of 𝑅⁺_𝑥, it is easy to verify that (see also [1]) 𝐸‖𝑠−𝑇𝑥‖²=tr𝑇−𝑇₀𝑅_𝑥𝑇−𝑇₀^𝑇+𝐸𝑠−𝑇₀𝑥². (3) Noting that the covariance matrix 𝑅_𝑥 is symmetric nonnegative definite, then it can be factorized as 𝑅_𝑥=𝑅_𝑥^1/2𝑅_𝑥^1/2𝑇. (4) Substituting (4) into (3) gives rise to 𝐸‖𝑠−𝑇𝑥‖²=tr𝑇−𝑇₀𝑅_𝑥^1/2𝑅_𝑥^1/2𝑇𝑇−𝑇₀^𝑇+𝐸𝑠−𝑇₀𝑥²=tr𝑇−𝑇₀𝑅_𝑥^1/2𝑇−𝑇₀𝑅_𝑥^1/2𝑇+𝐸𝑠−𝑇₀𝑥²=𝑇−𝑇₀𝑅_𝑥^1/22_𝐹+𝐸𝑠−𝑇₀𝑥²=𝑇𝑅_𝑥^1/2−𝑇₀𝑅_𝑥^1/22_𝐹+𝐸𝑠−𝑇₀𝑥²=𝑇₀𝑅_𝑥^1/2−𝑇𝑅_𝑥^1/22_𝐹+𝐸𝑠−𝑇₀𝑥², (5) since 𝐸{‖𝑠−𝑇₀𝑥‖²} is a constant, then min_{𝑇∈𝑅𝑚×𝑛(𝑇),rank=𝑑}𝐸‖𝑠−𝑇𝑥‖²=min_{𝑇∈𝑅𝑚×𝑛(𝑇),rank=𝑑}𝑇₀𝑅_𝑥^1/2−𝑇𝑅_𝑥^1/22_𝐹+𝐸𝑠−𝑇₀𝑥², (6) that is to say, minimizing 𝐸{‖𝑠−𝑇𝑥‖²} is equivalent to minimizing ‖𝑇₀𝑅_𝑥^1/2−𝑇𝑅_𝑥^1/2‖²_𝐹. Therefore, we can find the solution 𝑇^∗ of (1) by solving the minimization problem min_{𝑇∈𝑅𝑚×𝑛(𝑇),rank=𝑑}𝑇₀𝑅_𝑥^1/2−𝑇𝑅_𝑥^1/2_𝐹, (7) which can be summarized as the following low rank approximation problem:

Problem 1. Given two matrices 𝐴∈𝑅^𝑚×𝑛, 𝐵∈𝑅^𝑝×𝑛 and an integer 𝑑, 1≤𝑑<𝑚,𝑝, find a matrix 𝑋∈𝑅^𝑚×𝑝 of rank 𝑑 such that 𝐴−𝑋𝐵_𝐹=min_{𝑋∈𝑅𝑚×𝑝(𝑋),rank=𝑑}‖𝐴−𝑋𝐵‖_𝐹. (8)

In the last few years there has been a constantly increasing interest in developing the theory and numerical approaches for the low rank approximations of a matrix, due to their wide applications. A well-known method for the low rank approximation is the singular value decomposition (SVD) [5, 6]. When the desired rank is relatively low and the matrix is large and sparse, a complete SVD becomes too expensive. Some less expensive alternatives for numerical computation, for example, Lanczos bidiagonalization process [7], and the Monte Carlo algorithm [8] are available. To speed up the computation of SVD, random sampling has been employed in [9]. Recently, Ye [10] proposed the generalized low rank approximations of matrices (GLRAM) method. This method is proved to have less computational time than the traditional singular value decomposition-based methods in practical applications. Later, GLRAM method has been revisited and extended by Liu et al. [11] and Liang and Shi [12]. In some applications, we need to emphasize important parts and deemphasize unimportant parts of the data matrix, so the weighted low rank approximations were considered by many authors. Some numerical methods, such as Newton-like algorithm [13], left versus right representations method [14], and unconstrained optimization method [15], are proposed. Recently, by using the hierarchical identification principle [16] which regards the known matrix as the system parameter matrix to be identified, Ding et al. and Xie et al. present the gradient-based iterative algorithms [16–21] and least-squares-based iterative algorithm [22, 23] for solving matrix equations. The methods are innovational and computationally efficient numerical algorithms.

The common and practical method to tackle the low rank approximation Problem 1 is the singular value decomposition (SVD) (e.g. [1]). We briefly review SVD method as following. Minimizing (8) by a rank-𝑑 matrix 𝑋𝐵 is known [5, Page 69] to satisfy 𝑋𝐵=𝐴_𝑑=^𝑑_𝑖=1𝜎_𝑖𝑢_𝑖𝑣^𝑇_𝑖, (9) where 𝐴_𝑑 denotes rank-𝑑 singular value decomposition truncation, that is, if the following SVD holds 𝐴=^rank(𝐴)_𝑖=1𝜎_𝑖𝑢_𝑖𝑣^𝑇_𝑖, (10) then 𝐴_𝑑=∑^𝑑_𝑖=1𝜎_𝑖𝑢_𝑖𝑣^𝑇_𝑖. If the matrix 𝐵 is square and nonsingular, then by (9) we obtain that the solution of Problem 1 is 𝑋=𝐴_𝑑𝐵⁻¹=^𝑑_𝑖=1𝜎_𝑖𝑢_𝑖𝑣^𝑇_𝑖𝐵⁻¹. (11) The SVD method has two disadvantages as following: (1) it requires the matrix 𝐵 to be square and nonsingular; (2) in order to derive the solution (11), we must compute the inverse matrix of 𝐵, whose computation cost is very expensive.

In this paper, we develop a new method to solve the low rank approximation Problem 1, which can avoid the disadvantages of SVD method. We first transform Problem 1 into the fixed rank solution of a matrix equation and then use the generalized singular value decomposition (GSVD) to solve it. Based on these, we derive a sufficient condition for the existence of a solution of Problem 1, and the analytical expression of the solution is given. A numerical algorithm is proposed to compute the solution. Numerical examples are used to illustrate the numerical algorithm. The first one is artificial to show that the new algorithm is feasible to solve Problem 1, and the second is simulation, which shows that the new algorithm can be used to realize the image compression.

2. Main Results

In this section, we give a sufficient condition and an analytical expression for the solution of Problem 1 by transforming Problem 1 into the fixed rank solution of a matrix equation. Finally, we establish an algorithm for solving Problem 1.

Lemma 2. A matrix 𝑋∈𝑅^𝑚×𝑛 is a solution of Problem 1 if and only if it is a solution of the following matrix equation: 𝑋𝐵𝐵^𝑇=𝐴𝐵^𝑇,rank(𝑋)=𝑑. (12)

Proof. It is easy to verify that a matrix 𝑋∈𝑅^𝑚×𝑛 is a solution of Problem 1 if and only if 𝑋 satisfies the following two equalities simultaneously: 𝐴−𝑋𝐵_𝐹=min_{𝑋∈𝑅𝑚×𝑝}‖𝐴−𝑋𝐵‖_𝐹,rank(𝑋)=𝑑. (13)(14)
Since the normal equation of the least squares problem (13) is 𝑋𝐵𝐵^𝑇=𝐴𝐵^𝑇 (15) and noting that the least squares problem (13) and its normal equation (15) have the same solution sets, then (13) and (14) can be equivalently written as 𝑋𝐵𝐵^𝑇=𝐴𝐵^𝑇,rank(𝑋)=𝑑 (16) which also imply that Problem 1 is equivalent to (12).

Remark 3. From Lemma 2 it follows that Problem 1 is equivalent to (12), hence we can solve Problem 1 by finding a fixed rank solution of the matrix equation 𝑋𝐵𝐵^𝑇=𝐴𝐵^𝑇.

Now we will use generalized singular value decomposition (GSVD) to solve (12). Set 𝐶=𝐵𝐵^𝑇∈𝑅^𝑝×𝑝,𝐷=𝐴𝐵^𝑇∈𝑅^𝑚×𝑝. (17) The GSVD of the matrix pair (𝐶,𝐷) is given by (see [24]) 𝐶=𝑈Σ₁𝑊,𝐷=𝑉Σ₂𝑊, (18) where 𝑈∈𝑂^𝑝×𝑝, 𝑉∈𝑂^𝑚×𝑚, 𝑊∈𝑅^𝑝×𝑝 is a nonsingular matrix, 𝑘=rank([𝐶^𝑇,𝐷^𝑇]), 𝑟=rank(𝐶), 𝑡=rank(𝐶)+rank(𝐷)−rank([𝐶^𝑇,𝐷^𝑇]), and Σ₁=𝐼_𝐶000𝑆𝑟−𝑡^𝐶0𝑡00𝑂_𝐶000𝑡Σ𝑘−𝑟𝑝−𝑘𝑟−𝑡𝑝−𝑟,₂=𝑂_𝐷000𝑆𝑟−𝑡^𝐷0𝑡00𝐼_𝐷000𝑡𝑘−𝑟𝑝−𝑘𝑚−𝑘−𝑡+𝑟𝑘−𝑟 (19) are block matrices, with 𝐼_𝐶 and 𝐼_𝐷 are identity matrices, 𝑂_𝐶 and 𝑂_𝐷 are zero matrices: 𝑆_𝐶=diag𝛼₁,𝛼₂,…,𝛼_𝑡,1>𝛼₁≥𝛼₂≥⋯≥𝛼_𝑡𝑆>0,_𝐷=diag𝛽₁,𝛽₂,…,𝛽_𝑡,0<𝛽₁≤𝛽₂≤⋯≤𝛽_𝑡𝛼<1,_2𝑖+𝛽²_𝑖=1,𝑖=1,2,…,𝑡. (20)

By (17) and (18), we have 𝑋𝐵𝐵^𝑇=𝐴𝐵^𝑇⟺𝑋𝐶=𝐷⟺𝑋𝑈Σ₁𝑊=𝑉Σ₂𝑊⟺𝑋𝑈Σ₁𝑊−𝑉Σ₂𝑉𝑊=0⟺𝑉^𝑇𝑋𝑈Σ₁−Σ₂𝑊=0. (21) Set 𝑌=𝑉^𝑇𝑋𝑈, (22) and 𝑌 is partitioned as follows: 𝑌𝑌=¹¹𝑌₂₁𝑌₃₁𝑌𝑟−𝑡¹²𝑌₂₂𝑌₃₂𝑡𝑌₁₃𝑌₂₃𝑌₃₃𝑡𝑝−𝑟𝑚−𝑘−𝑡+𝑟𝑘−𝑟, (23) then 𝑉^𝑇𝑋𝑈Σ₁−Σ₂=𝑌Σ₁−Σ₂=𝑌₁₁𝑌₁₂𝑌₁₃𝑌₂₁𝑌₂₂𝑌₂₃𝑌₃₁𝑌₃₂𝑌₃₃𝐼_𝐶0000𝑆_𝐶0000𝑂_𝐶0−𝑂_𝐷0000𝑆_𝐷0000𝐼_𝐷0=𝑌₁₁𝑌₂₁𝑌₃₁𝑌𝑟−𝑡¹²𝑆_𝐶𝑌₂₂𝑆_𝐶−𝑆_𝐷𝑌₃₂𝑆_𝐶𝑡00−𝐼_𝐷000𝑡𝑘−𝑟𝑝−𝑘𝑚−𝑘−𝑡+𝑟𝑘−𝑟. (24)

Therefore, by (21) and (24), we have 𝑋𝐵𝐵^𝑇=𝐴𝐵^𝑇𝑌⟺𝑉¹¹𝑌₁₂𝑆_𝐶𝑌00₂₁𝑌₂₂𝑆_𝐶−𝑆_𝐷𝑌00₃₁𝑌₃₂𝑆_𝐶−𝐼_𝐷0𝑉𝑌𝑊=0⟺rank¹¹𝑌₁₂𝑆_𝐶𝑌00₂₁𝑌₂₂𝑆_𝐶−𝑆_𝐷𝑌00₃₁𝑌₃₂𝑆_𝐶−𝐼_𝐷0𝑊𝑌=0⟺rank¹¹𝑌₁₂𝑆_𝐶𝑌00₂₁𝑌₂₂𝑆_𝐶−𝑆_𝐷𝑌00₃₁𝑌₃₂𝑆_𝐶−𝐼_𝐷0=0⟺𝑘−𝑟=0,𝑌₁₁=𝑌₂₁=𝑌₃₁=𝑌₁₂=𝑌₃₂=0,𝑌₂₂=𝑆_𝐷𝑆_𝐶⁻¹, (25) that is to say, the matrix equation 𝑋𝐵𝐵^𝑇=𝐴𝐵^𝑇 has a solution if and only if 𝑘−𝑟=rank𝐶^𝑇,𝐷^𝑇−rank(𝐶)=rank𝐵𝐵^𝑇,𝐵𝐴^𝑇−rank𝐵𝐵^𝑇=0, (26) and according to (22), we know that the expression of the solution is 𝑋=𝑉𝑌𝑈^𝑇, (27) where 000𝑌=𝑆𝑟−𝑡^𝐷𝑆_𝐶⁻¹𝑡𝑌₁₃𝑌₂₃𝑌𝑝−𝑟𝑚−𝑡𝑡,₁₃∈𝑅^{(𝑚−𝑡)×(𝑝−𝑟)},𝑌₂₃∈𝑅^{𝑡×(𝑝−𝑟)}. (28)

By (26)–(28) and noting that 𝑌₁₃ and 𝑌₂₃ are arbitrary matrices, we have min_{𝑋∈𝑅𝑚×𝑝}rank(𝑋)=min_{𝑋∈𝑅𝑚×𝑝}rank𝑉𝑌𝑈^𝑇=min_{𝑌13∈𝑅(𝑚−𝑡)×(𝑝−𝑟),𝑌23∈𝑅𝑡×(𝑝−𝑟)}(𝑌)(𝐷)rank=𝑡=rank=rank𝐴𝐵^𝑇,max_{𝑋∈𝑅𝑚×𝑝}rank(𝑋)=max_{𝑌∈𝑅𝑚×𝑝}rank𝑉𝑌𝑈^𝑇=max_{𝑌13∈𝑅(𝑚−𝑡)×(𝑝−𝑟),𝑌23∈𝑅𝑡×(𝑝−𝑟)}(𝑌)(𝐶)(𝐷)𝐶rank=𝑡+min𝑝−𝑟,𝑚−𝑡=min𝑚,𝑝+𝑡−𝑟=min𝑚,𝑝+rank+rank−rank^𝑇,𝐷^𝑇(𝐶(𝐷)𝐶−rank)=min𝑚,𝑝+rank−rank^𝑇,𝐷^𝑇=min𝑚,𝑝+rank(𝐷)−rank(𝐶)=min𝑚,𝑝+rank𝐴𝐵^𝑇−rank𝐵𝐵^𝑇. (29)

Hence, if rank𝐵𝐵^𝑇,𝐵𝐴^𝑇−rank𝐵𝐵^𝑇=𝑘−𝑟=0,rank𝐴𝐵^𝑇≤𝑑≤min𝑚,𝑝+rank𝐴𝐵^𝑇−rank𝐵𝐵^𝑇, (30)(31) then (12) has a solution, and the expressions of the solution are given by (26)–(28), that is, 𝑋=𝑉𝑌𝑈^𝑇=𝑉00𝑌₁₃0𝑆_𝐷𝑆_𝐶⁻¹𝑌₂₃𝑈^𝑇, (32) where 𝑌₂₃∈𝑅^{𝑡×(𝑝−𝑟)} is an arbitrary matrix and 𝑌₁₃∈𝑅^{(𝑚−𝑡)×(𝑝−𝑟)} is chosen such that rank𝑌₁₃=𝑑−𝑡=𝑑−rank(𝐷)=𝑑−rank𝐴𝐵^𝑇. (33) And noting that the low rank approximation Problem 1 is equivalent to (12) (i.e. Lemma 2), then we obtain the following.

Theorem 4. If rank𝐵𝐵^𝑇,𝐵𝐴^𝑇−rank𝐵𝐵^𝑇=𝑘−𝑟=0,rank𝐴𝐵^𝑇≤𝑑≤min𝑚,𝑝+rank𝐴𝐵^𝑇−rank𝐵𝐵^𝑇, (34) then Problem 1 has a solution, and the expressions of the solution are given by 𝑋=𝑉𝑌𝑈^𝑇=𝑉00𝑌₁₃0𝑆_𝐷𝑆_𝐶⁻¹𝑌₂₃𝑈^𝑇, (35) where 𝑌₂₃∈𝑅^{𝑡×(𝑝−𝑟)} is an arbitrary matrix and 𝑌₁₃∈𝑅^{(𝑚−𝑡)×(𝑝−𝑟)} is chosen such that rank𝑌₁₃=𝑑−𝑡=𝑑−rank(𝐷)=𝑑−rank𝐴𝐵^𝑇. (36)

Remark 5. In contrast with (11), the solution expression (35) does not require the matrix 𝐵 to be square and nonsingular and does not need to compute the inverse of 𝐵.

Based on Theorem 4, we can establish an algorithm for finding the solution of Problem 1.

Algorithm 6. (1) Input the matrices 𝐴,𝐵 and the integer 𝑑; (2) make the GSVD of the matrix pair (𝐶,𝐷) according to (18); (3) choose 𝑌₂₃∈𝑅^{𝑡×(𝑝−𝑟)} and 𝑌₁₃∈𝑅^{(𝑚−𝑡)×(𝑝−𝑟)}, such that rank(𝑌₁₃)=𝑑−rank(𝐴𝐵^𝑇);(4) compute the solution 𝑋 according to (35).

3. Numerical Experiments

In this section, we first use a simple artificial example to illustrate that Algorithm 6 is feasible to solve Problem 1, then we use a simulation to show that Algorithm 6 can be used to realize the image compression. The experiments were done with MATLAB 7.6 on a 64-bit Intel Pentium Xeon 2.66 GHz with 𝑒_mach≈2.0×10⁻¹⁶.

Example 7. Consider Problem 1 with

,.𝐴=0.68840.58730.42360.42430.84831.04000.97780.15411.08690.82430.59980.59931.06951.30801.54600.35380.57090.39250.13440.34780.43180.53890.48660.10351.14660.89380.92780.52091.34101.50861.48490.25990.77520.58510.49720.39090.80550.93890.98320.1978𝐵=0.05540.63080.84680.94190.68910.28080.45580.77090.16550.47830.37020.75261.05480.43580.78350.33090.11500.38080.33160.59400.76770.31680.56450.29780.16930.50320.40010.79031.08900.44990.80730.35800.08690.42470.46080.64950.67790.27860.48290.41700.11420.57720.63510.88150.90400.37150.64210.5951 (37)

We make GSVD of the matrix pair (𝐶,𝐷)=(𝐵𝐵^𝑇,𝐴𝐵^𝑇) as follows: 𝐶=𝑈Σ₁𝑊,𝐷=𝑉Σ₂𝑊, (38) whereΣ𝑈=−0.75740.4215−0.05900.0130−0.2832−0.4060−0.1509−0.5893−0.3965−0.6437−0.1562−0.1845−0.1775−0.34610.8997−0.1454−0.0810−0.1071−0.1752−0.5847−0.16670.7506−0.1161−0.1510−0.3397−0.0817−0.0302−0.02510.9311−0.0970−0.4754−0.0814−0.0331−0.0218−0.09150.8703,𝑉=−0.28180.43220.09370.7418−0.4180−0.1814−0.6971−0.44170.0795−0.55030.9419−0.0182−0.03300.2445−0.22720.01660.5179−0.2279−0.5906−0.5751−0.0160−0.24240.8753−0.1863−0.3743,𝑊=−6.2034−5.7059−4.4070−5.9640−4.5789−6.1909−7.5279−8.6034−6.4574−8.9389−6.2317−8.37560.4175−0.13140.24530.4166−0.1876−0.73450.19280.66880.2243−0.6061−0.2199−0.2228−0.01010.2451−0.4669−0.02480.7439−0.40970.19260.3245−0.72130.2930−0.49100.1022₁,=0.9985000000.419200000000000000000000000,Σ₂=0000000000000000.0541000000.9079000. (39)(40)(41)(42)

It is easy to verify thatrank𝐵𝐵^𝑇,𝐵𝐴^𝑇−rank𝐵𝐵^𝑇=0,rank𝐴𝐵^𝑇=2,min𝑚,𝑝+rank𝐴𝐵^𝑇−rank𝐵𝐵^𝑇=5, (43) that is, if 2≤𝑑≤5, then the conditions of Theorem 4 are satisfied. Setting 𝑑=2∈[2,5], according to (35), we obtain that the solution of Problem 1 is 𝑋=𝑉0000000000000000000.054100.9985000000.9079𝑈0.41920000_𝑇,=−0.41210.52750.30620.52230.06030.0546−0.50570.70170.41170.69610.09590.0949−0.21750.28800.16800.28550.03570.0338−0.50080.73890.43670.73390.11260.1166−0.33410.47930.28230.47580.06960.0708. (44)(45)

Setting 𝑑=4∈[2,5], according to (35), we obtain that the solution of Problem 1 is 𝑋=𝑉0010000001000000000.054100.9985000000.9079,𝑈0.41920000_𝑇=.−0.38980.3610−0.01020.89370.05790.0549−0.50401.22240.34980.20320.11880.1155−0.2733−0.07361.01810.11480.00770.0030−0.49500.39890.37641.11980.09910.1052−0.33630.64160.30320.29650.07620.0764 (46)

Example 7 shows that Algorithm 6 is feasible to solve Problem 1. However, the SVD method in [1] cannot be used to solve Example 7, because 𝐵 is not a square matrix.

Example 8. We will use the generalized Karhunen-Loeve transform, based on Algorithm 6 and SVD method in [1], respectively, to realize the image compression. Figure 1(a) (see page 3) is the test image which has 256×256 pixels and 256 levels on each pixel. We separate it into 32×32 blocks such that each block has 8×8 pixels. Let 𝑓^(𝑘,𝑙)_𝑖,𝑗 and 𝑛^(𝑘,𝑙)_𝑖,𝑗 (𝑖,𝑗=0,1,2,…,7;𝑘,𝑙=0,1,2,…,31) be the values of the image and a Gaussian noise (generated by Matlab function 𝑖𝑚𝑛𝑜𝑖𝑠𝑒) at the (𝑖,𝑗)th pixel in the (𝑘,𝑙)th block, respectively. For convenience, let 𝑎=𝑖+8𝑗, 𝑝=𝑘+32𝑙, and the (𝑖,𝑗)th pixel in the (𝑘,𝑙)th block be expressed as the 𝑎th pixel in the 𝑝th block (𝑎=0,1,2,…,63;𝑝=0,1,…,1023). We can also express 𝑓^(𝑘,𝑙)_𝑖,𝑗 and 𝑛^(𝑘,𝑙)_𝑖,𝑗 as 𝑓_𝑎^(𝑝) and 𝑛_𝑎^(𝑝), respectively.
The test image is processed on each block. Therefore, we can assume that the blocked image space is 64-𝐷 real vector space 𝑅⁶⁴. The 𝑝th block of the original image is expressed by the 𝑝th vector: 𝑠^𝑝=𝑠^𝑝₀,𝑠^𝑝₁,…,𝑠^𝑝₆₃^𝑇. (47) Hence the original image is expressed by 1024 64-𝐷 vectors {𝑠^𝑝}¹⁰²³_𝑝=0. The noise is similarly expressed by {𝑛^𝑝}¹⁰²³_𝑝=0, where 𝑛^𝑝=𝑛^𝑝₀,𝑛^𝑝₁,…,𝑛^𝑝₆₃^𝑇. (48) Figure 1(b) is the noisy image {𝑥^𝑝}¹⁰²³_𝑝=0, where 𝑥^𝑝=𝑠^𝑝+𝑛^𝑝,𝑝=0,1,…,1023. (49) By (47), (49), (2), (4) and the definition of covariance matrix, we get 𝑇₀ and 𝑅_𝑥^1/2 of (7). Then we use Algorithm 6 and SVD method in [1] to realize the image compression respectively, and the experiment results are in pages 4 and 5.
Figure 2 illustrates that Algorithm 6 can be used to realize image compression. Although it is difficult to see the difference between Figures 2 and 3, which are compressed by SVD method in [1], from Table 1 we can see that the execution time of Algorithm 6 is less than that of SVD method at the same rank. This shows that our algorithm outperforms the SVD method in execution time.

Table 1 

Execution time for deriving Figures 2(a)–3(c).

(a)

(b)

Figure 1 

(a) Original image; (b) noisy image.

(a)

(b)

(c)

Figure 2 

Image compression by Algorithm 6 with different rank 𝑑: (a) 𝑑=40; (b) 𝑑=30; and (c) 𝑑=20.

(a)

(b)

(c)

Figure 3 

Image compression by SVD method with different rank 𝑑: (a) 𝑑=40; (b) 𝑑=30; and (c) 𝑑=20.

4. Conclusion

The low rank approximation Problem 1 arising in the generalized Karhunen-Loeve transform is studied in this paper. We first transform Problem 1 into the fixed rank solution of a matrix equation and then use the generalized singular value decomposition (GSVD) to solve it. Based on these, we derive a sufficient condition for the existence of a solution, and the analytical expression of the solution is also given. Finally, we use numerical experiments to show that new algorithm is feasible and effective.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (11101100; 11226323; 11261014; and 11171205), the Natural Science Foundation of Guangxi Province (2012GXNSFBA053006; 2013GXNSFBA019009; and 2011GXNSFA018138), the Key Project of Scientific Research Innovation Foundation of Shanghai Municipal Education Commission (13ZZ080), the Natural Science Foundation of Shanghai (11ZR1412500), the Ph.D. Programs Foundation of Ministry of Education of China (20093108110001), the Discipline Project at the corresponding level of Shanghai (A. 13-0101-12-005), and Shanghai Leading Academic Discipline Project (J50101).