CN108664448A - A kind of higher-dimension sparse vector reconstructing method based on IQR_OMP - Google Patents

Abstract

The invention discloses a kind of higher-dimension sparse vector reconstructing method based on IQR_OMP.The characteristics of will calculating least square solution is often walked according to OMP algorithms, increment IGR_OMP algorithms is proposed based on Greville recursive processes, and obtain a series of useful recursion properties, is effectively improved calculating speed；On the basis of analyzing IGR_OMP algorithms, the relationship decomposed using Greville recursive algorithms and QR establishes the increment IQR_OMP algorithms on QR decomposition bases, which effectively reduces amount of calculation, improve the validity of algorithm solution large-scale problem；Higher-dimension sparse vector reconstructs simply exactly provides the most sparse coefficient expression of high dimension vector on excessively complete system of vectors.

Description

A kind of higher-dimension sparse vector reconstructing method based on IQR_OMP

Technical field

The present invention relates to a kind of higher-dimension sparse vector reconstructing method based on IQR_OMP.

Background technology

Compressed sensing (compressive Sensing, abbreviation CS) is that one of image or field of signal processing is new Key technology.Its essence be using the information far below signal vector or image vector dimension, effectively determine this higher-dimension to Most sparse coefficient of the amount on excessively complete dictionary indicate, key problem is just to solve for the sparse of sub- alignment sexual system Solution, with reconstruction signal or image vector (bibliography：Yan Jingwen etc., compressed sensing and application, 2015, National Defense Industry Press, Beijing .E.J.Candes, M.Wakin, " people hearing without listening " An introduction to compressive sampling,IEEE Signal Processing Magazine,2008,25(2),21-30.)。

Higher-dimension sparse vector reconstruct be signal processing, image processing field a particularly significant problem.Classical OMP (Orthogonal Matching Pursuit, orthogonal matching pursuit) is although complexity of the algorithm relative to other non-greedy algorithms It spends relatively low, but still has a distance apart from real-time signal processing in practical applications.

Invention content

The present invention in view of the deficiencies of the prior art, discloses a kind of based on IGR_OMP (Incremental Greville Recursion based Orthogonal Matching Pursuit, orthogonal matching of the incremental based on Greville recursion Tracing algorithm) higher-dimension sparse vector reconstructing method, include the following steps：

If matrix A=[a₁, a₂..., a_K]∈R^M×K, K ＜ ＜ M, rank (A)=K (R^M×K：M × K ties up real number matrix set, a_i(i=1,2 ..., K)：I-th of real number column vector in dictionary matrix A, rank (A)：The order of matrix A)；

Remember A₁=[a₁], A_k=[A_k-1,a_k], k=2,3 ..., K (a₁：First column vector of matrix A, A₁：1st iteration Matrix A, A_k：There are k-th of matrix A of k column vector, A_k-1：- 1 Iterative Matrix A of kth)；

A is calculated by Greville recursive algorithms⁺, i.e. generalized inverse (the bibliography of matrix A：Ben-Isreal A., Greville T.N.E,Generalized Inverse Theory and Applications,Wiley,New York, 1974.)：

Preliminary examination value

Recursion is walked to k=2,3 ..., K

b_k=a_k-A_k-1d_k；

Wherein d_kIt is to solve forThe vector replaced temporarily, b_kIt is dimension and a_kIdentical solutionThe column vector replaced temporarily. The characteristics of this algorithm isByRecursion obtains, that is, fromGradually calculateAnd It is not all directly to calculate each iterationIt is for the sparse recovering step of increment IGR_OMP high dimension vectors in compressed sensing：

Step A1, input observation data vector y ∈ R^MWith dictionary matrix A=[a₁, a₂..., a_N]∈R^M×N, enable initial support Collect indexEnable initial residual vector r₀Equal to observation data vector y, i.e. r₀=y, enables k=1.Wherein R^MIt is tieed up for M × 1 real Ordered series of numbers vector set, R^M×NFor M × N-dimensional real number matrix set, a_iRepresent i-th of real number column vector in dictionary matrix A, i=1, 2…,N；

Step A2, ask in dictionary matrix A with -1 iteration residual vector r of kth_k-1Kth time iteration strongest correlation row j_k：

j_k∈arg max_j|(r_k-1, a_j) |, Ω_k=Ω_k-1∪{j_k,

Wherein Ω_kFor kth time iteration supported collection index, Ω_k-1For -1 iteration supported collection index of kth, a_jFor dictionary matrix A In j-th of real number column vector；

Step A3, by -1 iteration support index Ω of kth_k-1Corresponding dictionary group inverse matrices matrix Incremental computations kth time iteration support index Ω_kCorresponding dictionary group inverse matrices matrix

If k=1, j_kFor j₁, then vectorialGeneralized inverse be equal toGeneralized inverse, i.e., Turn to step A4；

If k ≠ 1, sequence executes following calculate：

C=b^T/(b^TB),

WhereinIt is dictionary matrixMiddle jth_kA real number column vector,It is vectorGeneralized inverse matrix, d is to face When a vector replacing, b is the dimension replaced temporarily and vectorThe dimensional vectors of identical M × 1, c are the 1 × M replaced temporarily Dimensional vector.

Step A4 solves minimization problemSolution x_k：

WhereinRepresentation vectorEuclid norm；

Step A5, update kth time iteration residual vector

Step A6 judges whether to meet stopping criterion, that is, judge kth time iteration residual vector Euclid norm whether Meet ‖ r_k‖₂(ε is a very small number to≤ε, such as：10^-6), if it is satisfied, executing step A7, k is otherwise updated to k+1, And step A2 is repeated to step A5；

Step A7 exports sparse coefficient vector x ∈ R^N：

Wherein x (i) is i-th of element in sparse coefficient vector x, x_k(i) it is vector x_kIn i-th of element, R^NIt is N × 1 ties up real vector.

b_k=0 necessary and sufficient condition is a_k∈R(A_k-1)=span { a₁, a₂..., a_k-1,

That is a_kWith a₁, a₂..., a_k-1Linear correlation, A_k-1=[a₁, a₂..., a_k-1](a_kRepresenting matrix A_kIn arrange for k-th to Amount, R (A_k-1) represent with matrix A_k-1K-1 column vector generate space, span { a₁, a₂..., a_k-1Represent with vectorial a₁, a₂..., a_k-1The vector space of generation, A_k-1For -1 iteration A matrix of kth,For A_k-1Generalized inverse matrix)；

b_k∈R(A_k-1)^⊥It is a_kIn R (A_k-1) rectangular projection in orthocomplement, orthogonal complement, if b_k=0, indicate a_k∈R(A_k-1), then a_k With R (A_k-1) linear correlation (I_MIt is the unit matrix of M × M,Representing matrix A_k-1Generalized inverse matrix, R (A_k-1)^⊥Indicate life At space R (A_k-1) the orthogonal complement space).

Indicate b_kIt is a_kIn R (A_k-1)^⊥On rectangular projection, then have (a_i(i=1, 2 ..., k-1) it is matrix A_k-1In the i-th column vector)：

(b_k, a_i)=0, i=1,2 ..., k-1 (4.7)

That is b_kWith a₁, a₂..., a_k-1It is orthogonal.

The invention also discloses one kind being based on IQR_OMP (Incremental QR decomposition based Orthogonal Matching Pursuit, the orthogonal matching pursuit algorithm that incremental is decomposed based on QR) higher-dimension sparse vector Reconstructing method includes the following steps：

Step B1, input observation data vector y ∈ R^mWith dictionary matrix A=[a₁, a₂..., a_n]∈R^m×n；Enable initial support Collect indexInitial residual vector is enabled to be equal to observation data vector, i.e. r₀=y enables k=1, wherein R^mIndicate that the dimensions of m × 1 are real Ordered series of numbers vector, R^m×nIndicate that m × n ties up real number matrix, a_iIt represents i-th of m × 1 in dictionary matrix A and ties up real number column vector, i=1, 2…,n；

Step B2, ask in dictionary matrix A with -1 iteration residual vector r of kth_k-1Kth time iteration strongest correlation row j_k：

j_k∈arg max_j|(r_k-1, a_j) |, Ω_k=Ω_k-1∪{j_k,

Wherein Ω_kFor kth time iteration supported collection index, Ω_k-1For -1 iteration supported collection index of kth, a_jFor dictionary matrix A In jth row real number column vector；

Step B3, by the Q of -1 iteration of kth_k-1And R_k-1The Q of recurrence calculation kth time iteration_kAnd R_k：

If k=1, j_k=j₁, thenQ₁The first of matrix is classified as Q₁(:, 1) or q₁It is equal toI.e.Step B4 is turned to,

If k ≠ 1, following calculate is executed successively：

q_k=q₁/α_k,

Q_k=[Q_k-1:q_k],It is wherein vectorialIt is the jth of dictionary matrix A_kA column vector, q_kIt is K-th of column vector of matrix Q, β_kInterim by kth time iteration replaces vector, α_kInterim by kth time iteration replaces numerical value；

Step B4, update kth time iteration residual vector

Wherein q_kIt is kth time Iterative Matrix Q_kK-th of column vector；

Step B5 judges whether to meet stopping criterion, that is, judges whether the Euclid norm of kth time residual vector is full Sufficient ‖ r_k‖₂(ε is a very small number to≤ε, such as：10^-6), if it is satisfied, executing step B6；Otherwise k is updated to k+1, and Step B2 is repeated to step B4；

Step B6 solves minimization problemSolution x_k：

WhereinFor kth time iteration support index Ω_kCorresponding dictionary matrix, j₁, j₂..., j_k∈Ω_k,Representation vectorEuclid norm；It is matrix R_kGeneralized inverse,It is matrix Q_kTransposition；

Step B7 exports sparse coefficient vector x：

Wherein x (i) indicates i-th of element in sparse coefficient vector x, x_k(i) vector x is indicated_kIn i-th of element.

It is obtained using Greville recurrence Relations：

Then have：Wherein d_kIt is to solve forIt replaces temporarily Vector, b_kIt is dimension and a_kIdentical solutionThe column vector replaced temporarily.

(QR Decomposition) is decomposed according to QR, if having calculated the matrix A of -1 execution of kth_k-1：

A_k-1=Q_k-1R_k-1 (4.9)

Wherein A_k-1Kth -1 Iterative Matrix A, Q_k-1=[q₁,q₂,…q_k-1]∈R^M×(k-1)Row directly hand over matrix (q_i(i= 1,2 ..., k-1) it is matrix Q_k-1The i-th column vector), q_k-1It is matrix Q_k-1- 1 column vector R of kth_k-1For (k-1) × (K-1) Upper triangular matrix, consider kth step calculating：

Wherein β_k∈R^(k-1)×1, α_kFor a scale, (I is pushed away：Unit matrix)：

Due toThusWhereinFor matrix R_kInverse square Battle array,

So：

Equally often step calculates residual vector r_k=s-A_kx_kIt needs to calculate a matrix and multiplication of vectors (s is reconstruct vector), And r_kOften step needs to be used in selection and r in A_kInner product adds to A by the column vector of maximum absolute value_kMiddle formation A_k+1, at this moment for Recursion QR has for decomposing：r_kIt is y in R (A_k)^⊥On rectangular projection, i.e.,：

Advantageous effect：

Two increment OMP methods of IGR_OMP methods and IQR_OMP of the sparse reconstruct of high-dimensional signal vector given herein There are following features：

1. calculating speed is fast.For matrix A_k=[a₁, a₂..., a_k]∈R^M×k, y ∈ R^M, rank (A_k)=k is corresponding min_x‖y-A_kThe solution of x ‖ isIf decomposing to obtain by QRIt is traditional for each k OMP methods, calculation amount are aboutSecondary flop (bibliography：L.N.Trefethen,D.Bau, Numerical Linear Algebra,Chinese edition copyright 2006by Posts and Telecom Press).And delta algorithm IQR_OMP algorithms are aboutTherefore the workload of incremental algorithm is substantially better than tradition OMP algorithms, especially when higher-dimension problem.

2. from the point of view of computational stability, due toAnd incremental IQR_OMP Method design is cond₂(A_k), rather than its square.Therefore work as A_kConditional number it is slightly larger when, normal equation algorithm is not only It is computationally intensive, it can also seriously affect the stability of algorithm.Therefore the present invention is substantially better than in computational accuracy and in calculating speed Normal equation method.

3. the notional result about the sparse reconstruct of signal vector is difficult to be examined during actually calculating at present. And the r in the delta algorithm IQR_OMP and IGR_OMP provided in text_k,kAnd b_k, can be used as and examine a_kA relative to front₁, a₂,…,a_k-1Measurement (the bibliography of linear independent degree：Zhao Jinxi, the Huang methods of compatible linear equation group and its pushes away Extensively, institution of higher education calculates mathematics journal, 01 phase, pp 8-17,1981/4/2.).It realizes and examines vector during calculating Linear independent degree between system, this is very important.The characteristics of incremental algorithm that the present invention provides is to substantially reduce meter Expense is calculated, algorithm has very strong robustness, and currently selected atom can be provided during calculating relative to front The linear independence degree of atom.Numerical result shows that algorithm given herein has apparent advantage.

Description of the drawings

The present invention is done with reference to the accompanying drawings and detailed description and is further illustrated, it is of the invention above-mentioned or Otherwise advantage will become apparent.

Fig. 1 is IGR_OMP algorithm flow charts.

Fig. 2 is IQR_OMP algorithm flow charts.

Specific implementation mode

The present invention will be further described with reference to the accompanying drawings and embodiments.

In the sparse restructing algorithm of the higher-dimension signal vector of compressed sensing, greedy (Greedy) algorithm is with its lower calculating Complexity and faster convergence rate are paid attention to, and orthogonal matching pursuit algorithm (OMP) is one algorithm of most important one (bibliography：Joel A.Tropp and Anna C.Gilbert.Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit[J].IEEE Transactions on Information Theory,VOL.53,NO.12,DECEMBER 2007.)。

Classical orthogonal matching pursuit algorithm (OMP, Orthogonal Matching Pursuit) is as follows：

Input：Observe data vector y ∈ R^M(R^M：M ties up real number column vector set) and dictionary matrix A=[a₁, a₂..., a_N] ∈R^M×N(a_i(i=1,2 ..., N) is the i-th column vector of matrix A, R^M×N：M × N-dimensional real number matrix).

Output：Sparse signal vector x ∈ R^N(R^N：N-dimensional real number column vector set).

Step1. it initializes：Enable initial support collection indexPreliminary examination residual vector is enabled to be equal to observation data vector, i.e., r₀=y, enables k=1.

Step2. it recognizes：Ask in matrix A with -1 iteration residual vector r of kth_k-1Kth time iteration strongest correlation row j_k：

Wherein Ω_kFor kth time iteration supported collection index, Ω_k-1For -1 iteration supported collection index of kth, a_j(j=1,2 ... K) j-th of column vector in matrix A.

Step3. it calculatesSolution

WhereinIt is support index Ω_kCorresponding dictionary matrix,a_i(i=j₁, j₂... j_k∈Ω_k) it is matrixI-th column vector.

Step4. update kth time iteration residual vector

If Step5. some stopping criterion meets, i.e. the Euclid norm of kth time iteration residual vector meets ‖ r_k‖₂≤ε (ε is a very small number, such as：10^-6), then stop iteration and turns Step6；

Otherwise k ← k+1 is enabled, and repeats Step2 to Step4.

Step6. output factor vector x (x (i)：Indicate i-th of element of output factor vector x, x_k(i)：Indicate vector x_k I-th of element)：

Due to non-zero entry number K≤M ＜ ＜ N of solution vector, step3 is actually the overdetermination for solving a M × k in the algorithm Linear least squares minimization problem.That is it is to seek the sparse solution of sub- alignment sexual system, and in fact OMP algorithms are on the whole It is realized by gradually acquiring the determined linear least square problem of M × k (k=1,2 ..., K).That is it is calculated with greedy Method seeks the sparse reconstruction of higher-dimension signal vector, is exactly the column vector of A corresponding to the nonzero element asked by x in essence The least square solution for the determined linear system that (atom) is formed, it is understood that this point is critically important.There are two crucial for OMP algorithms Step：

1. for matrix (dictionary) A ∈ R^M×NDetermine row (atom) index to match with nonzero element in x.According to greediness The principle of algorithm in kth step is found out and r in the N row of A_k-1=y-Ax_k-1Vector (the bibliography of inner product maximum absolute value： S.Mallat and Z.Zhang,Matching pursuits with time-frequency dictionaries,IEEE Trans.Signal Process, 1993,41, pp.3397-3415.), new supported collection is formed, here it is step2；

2. for large-scale problem, the work needed for least square solution is solved with common normal equation method in OMP algorithms Amount will often account for nearly the 80% of entire workload.Therefore finding the more efficiently algorithm for solving the least square solution just becomes A key problem in OMP algorithms.

Noted above, the sparse reconstruct of signal vector is to determine that a large-scale high-order Asia is fixed for its problem in science The most sparse solution of linear system.And OMP algorithms are actually introduced supported collection by column and are gradually solved such as next determined linear The least square solution of system：

And common OMP algorithms are normal equation methods, are calculated when kth walks iteration：

For solving the least square solution (3.3) of overdetermined linear system, due to iterations K ＜ ＜ M ＜ ＜ N. See that it is suitable effective for other direct methods to obtain (3.4) by normal equation method in this way.But here There are one basic problem, the coefficient matrix in (3.3)It is that each iteration increases by a row, and iteration will solve every time (3.3), the calculation amount accumulated is just very big.To general middle-size and small-size problem, also unobvious are influenced, but to large-scale high-order Problem, here it is a prodigious problems.How the information of OMP algorithm successive iteration is effectively utilized, this is construction of the present invention The starting point of incremental new algorithm.

In order to more intuitively describe the problem, if in the present invention determined linear problem (3.3) to be solvedNamely 2000 step of iteration, then the accumulation of usage equation method needs to calculate 1667.6 seconds, and present invention profit It is calculated and has only been used 137.36 seconds with the accumulation of the new algorithm of the information structuring of previous iteration, same problem speed improves 12 times. In the 2000th iteration, normal equation method has been used 4.5273 seconds, and the method for the present invention has only been used 0.2427 second.

Front, which is talked about, to discuss the existence of Large Scale Sparse solution using the concept of spark or the cross correlation value of A.But It has to be observed that these conditions are all less convenient for examining.In fact due to the present invention solve be one M × k determined linear system, institute Just to have：

Theorem 3.1 sets sub- alignment sexual system Ax=y, wherein A ∈ R^M×N, x ∈ R^N, y ∈ R^M, rank (A)=M, ‖ x ‖₀=K and K ＜ ＜ M ＜ ＜ N, if A=[a₁, a₂..., a_N] arbitrary K row composition system of vectors all linear independences, then certainly exist satisfaction The solution vector x of sparse condition so that ‖ x ‖₀=K.

This condition seems to be difficult verification, but since the coefficient matrix of this determined linear system is introduced by column, As long as ensureing that the system of vectors for the vector composition that this is introduced by column is linear independence in fact.How during calculating Can verify that newly into vector with into this point of the system of vectors linear independence of supported collection with regard to particularly important, and traditional OMP Algorithm does not accomplish this point.The characteristics of how being introduced by column using this considers and algorithm for design that this is that the present invention considers new The motivation of algorithm.

As depicted in figs. 1 and 2, the invention discloses a kind of higher-dimension sparse vector reconstructing method and one based on IGR_OMP Higher-dimension sparse vector reconstructing method of the kind based on IQR_OMP.

The principle of increment OMP methods：

OMP algorithms are substantially one determined linear systems of solution under certain condition, and the spy of this determined linear system Point is that gradually the selection maximum column vector of energy enters supported collection.Namely at the kth iteration, selection energy it is maximum arrange to M × k overdetermined systems of amount and the original k-1 Column vector groups Cheng Xin selectedWhereiny∈R^M, so After seek its least square solution.The essence that just can more catch OMP methods is understood and analyzed with this viewpoint, can more find out tradition It is insufficient existing for OMP methods.The each iteration of OMP algorithms will calculated sub-matrix from the beginning generalized inverse.Either directly adjust It is still directly calculated using formula with canonical functionMany additional calculating will be increased, especially to large-scale problem with The increase of iterations, this problem are just all the more prominent.And the characteristics of delta algorithm kth time iteration has been obtained using preceding an iteration To information, to greatly improve the validity of algorithm.Two increment OMP algorithms are given below.

IGR_OMP algorithms：

If A=[a₁, a₂..., a_K]∈R^M×K, K ＜ ＜ M, rank (A)=K (wherein a_i(i=1,2 ..., K) it is matrix A I-th of column vector, R^M×K：M × K ties up real number matrix set)；

Remember A₁=[a₁], A_k=[A_k-1,a₁], k=2,3 ..., K (a₁：First column vector of matrix A, A₁：1st iteration Matrix A, A_k：Kth time Iterative Matrix A, A_k-1：- 1 Iterative Matrix A of kth)；

A is calculated by Greville recursive algorithms⁺(bibliography：Ben-Isreal A.,Greville T.N.E, Generalized Inverse Theory and Applications,Wiley,New York,1974.)：

Preliminary examination value

Recursion is walked to k=2,3 ..., K

b_k=a_k-A_k-1d_k；

The characteristics of this algorithm isByRecursion obtains, wherein d_kIt is to solve forThe vector replaced temporarily, b_kIt is dimension Degree and a_kIdentical solutionThe column vector replaced temporarily.Therefore there are following IGR_OMP algorithms：

Step A1, input observation data vector y ∈ R^M(R^M：M × 1 ties up real number column vector) and dictionary matrix A=[a₁, a₂..., a_N]∈R^M×N(a_i(i=1,2 ..., N) is the i-th column vector of matrix A, R^M×N：M × N-dimensional real number matrix), enable initial branch Support collection indexInitial residual vector is enabled to be equal to input observation data vector, i.e. r₀=y, enables k=1；

Step A2, ask in dictionary matrix A with -1 iteration residual vector r of kth_k-1Kth time iteration strongest correlation row j_k (Ω_k：Kth time iteration supported collection index, Ω_k-1：- 1 iteration supported collection index of kth, a_j：J-th of column vector in matrix A)：

Step A3, by -1 iteration supported collection index Ω of kth_k-1Corresponding matrixGeneralized inverseIncrement meter Calculate kth time iteration supported collection index Ω_kCorresponding matrixGeneralized inverse(Represent matrixJth_kRow)：

If k=1, j_k=j₁, then vectorialGeneralized inverse be equal to vectorGeneralized inverse, i.e., Step A4 is turned to,

If k ≠ 1 executes the following steps successively：

C=b^T/(b^TB),

Wherein d be it is interim replace vector, b be it is interim replace withThe identical column vector of dimension, c be the row that replaces temporarily to Amount.

Step A4, minimization problemSolution

Wherein

Step A5, update kth time iteration residual vector

Step A6, judges whether to meet stopping criterion, i.e. the Euclid norm of kth time iteration residual vector meets ‖ r_k‖₂ (ε is a very small number to≤ε, such as：10^-6), if it is satisfied, executing step A7, k is otherwise updated to k+1, and repeat step A2 to step A5；

Step A7 exports sparse coefficient vector x ∈ R^N(R^N：N × 1 ties up real number column vector, x (i)：Export sparse coefficient I-th of element of vector x, x_k(i)：Step A4 minimization problem solution vectors x_kI-th of element)：

Arrange the property of Greville recurrence methods：

Property 4.1b_k=0 necessary and sufficient condition is a_k∈R(A_k-1)=span { α₁, a₂..., a_k-1, that is, a_kWith a₁, a₂..., a_k-1It is linearly related.

This is because：

This illustrates b_k∈R(A_k-1)^⊥, it is a_kIn R (A_k-1) rectangular projection in orthocomplement, orthogonal complement, if b_k=0, it is meant that a_k ∈R(A_k-1), therefore a_kWith R (A_k-1) linearly related, otherwise also set up.

The ‖ b in actually calculating_k‖₂It is a critically important index, it reflects a_kWith R (A_k-1) linearly related degree.

In compressed sensing problem, the hypothesis by perception matrix is not in b_k=0 the case where.But ‖ b_k‖₂This index is non- Often important, it reflects a_kWith R (A_k-1) linearly related degree.If ‖ b_k‖₂≠ 0, but very little just illustrates a_kWith a₁, a₂..., a_k-1Very weak close to linear correlation, that is, linear independence degree, here it is the unstable situations of so-called numerical value.Its Real this point is extremely important, it solves one during calculating to examine a_kWith the front system of vectors a₁, a₂..., a_k-1 Linear independence degree, this is the big problem in a numerical computations.

Property 4.2That is b_kIt is a_kIn R (A_k-1)^⊥On rectangular projection, therefore have：

(b_k, a_i)=0, i=1,2 ..., k-1 (4.7)

That is b_kWith a₁, a₂..., a_k-1It is orthogonal.

IQR_OMP algorithms：

If having calculated：

A_k-1=Q_k-1R_{k_1} (4.9)

Wherein Q_{k_1}=[q₁,q₂,…q_k-1]∈R^M×(k-1)Row directly hand over matrix, R_k-1For upper three angular moment of (k-1) × (k-1) Battle array.The calculating of kth step is considered below：

Wherein β_k∈R^(k-1)×1, α_kFor a scale.It pushes away：

Due toThus

So：

Equally often step calculates residual r_k=s-A_kx_kIt needs to calculate a matrix and multiplication of vectors (s：Need to reconstruct to Amount), and r_kOften step needs to be used in selection and r in A_kThe maximum column vector of inner product adds to A_kMiddle formation A_k+1.At this moment for recursion QR has for decomposing

Property 4.4 is in QR_OMP algorithms, r_kIt is y in fact in R (A_k)^⊥On rectangular projection, i.e.,：

K=1,2 ..., K

Residual vector r is calculated in (4.9) formula of utilization_kIt brings great convenience, avoids matrix and be directly multiplied with vector, often Step calculate inner product can, due to avoiding calculating intermediate quantity { x in recycling inside_k, this can reduce large-scale problem Amount of calculation.Therefore have：

IQR-OMP algorithms：

Step B1, input observation data vector y ∈ R^m(R^m：M × 1 ties up real number column vector set) dictionary matrix A=[a₁, a₂..., a_n]∈R^m×n(a_i(i=1,2 ... n) be matrix A the i-th column vector, R^m×n：M × n ties up real number matrix set)；It enables just Beginning supported collection indexInitial residual vector is enabled to be equal to observation data vector, i.e. r₀=y, enables k=1；

Step B2, ask in dictionary matrix A with -1 iteration residual vector r of kth_k-1Strongest correlation row j_k：

j_k∈arg max_j|(r_k-1, a_j) |, Ω_k=Ω_k-1∪{j_k}；

Wherein Ω_kIt is kth time iteration supported collection index, Ω_k-1It is -1 iteration supported collection index of kth, a_jIt is in matrix A J-th of column vector.

Step B3, by -1 iteration Q of kth_k-1And R_k-1Recurrence calculation kth time iteration Q_kAnd R_k：

If k=1, j_k=j₁, thenMatrix Q₁First be classified as Q₁(:, 1) or q₁It is equal toI.e.

Step B4 is turned to,

If k ≠ 1 executes following formula successively：

q₁=a_jk-Q_k-1β_k,

q_k=q₁/α_k,

Q_k=[Q_k-1:q_k],WhereinIt is jth in matrix A_kA column vector, q_kIt is the of matrix Q K column vector, β_k、α_kThe interim substitute element executed for kth time；

Step B4, update kth time iteration residual vector

Step B5, judges whether to meet stopping criterion, i.e. the Euclid norm of kth time iteration residual vector meets ‖ r_k‖₂ (ε is a very small number to≤ε, such as：10^-6), if it is satisfied, executing step B6；

Otherwise k is updated to k+1, and repeats step B2 to step B4；

Step B6, minimization problemSolution：

It providesIt is support index Ω_kCorresponding dictionary matrix,a_i(i=j₁, j₂... j_k∈Ω_k) it is matrixI-th column vector.

Step B7 exports sparse coefficient vector x (x (i)：Export i-th of element of sparse coefficient vector x, x_k(i)： Step B6 acquires solution x_kI-th of element)：

It can be obtained using Greville recurrence Relations：

Therefore have：

The process of theorem 4.5Grville recurrence calculation M-P generalized inverses is consistent with QR recursion.In fact in theorem 3 Under the conditions of, have

Embodiment

The computing environment of this paper design values experiment is DELL-PC WORK GROUP Intel (R) Core (TM) i7- 4790 3.6Ghz, memory 8.0GB, Windows7 Ultimate, software for calculation are MATLAB R2010b.

A ∈ R in Case1^2000×5000, K=600；

A ∈ R in Case2^3000×6000, K=1000；

A ∈ R in Case3^10000×15000, K=3000；

The calculating time of increment OMP algorithms 1-D is compared as follows shown in table 1：

Table 1

Situation type	IQR_OMP	IGR_OMP	OMP
				Case1	7.6331	9.112	30.8978
Case2	29.2647	34.8170	164.7083
				Case3	810.4921	962.1627	10386.0

Three above algorithmic descriptions are as follows：

Given herein establish is decomposed and algorithm on the basis of new residual recurrence formula in QR when 1.IQR_OMP algorithms；

It is given herein when 2.IGR_OMP algorithms to establish in Greville recurrence formula and new residual recurrence formula basis On OMP algorithms；

3.OMP algorithms are exactly the OMP algorithms of classics, the wherein form of generalized inverse normal equation.I.e.Simultaneously it is to be appreciated that the above same type, the initial data A ∈ R used in algorithms of different^M×NWith y ∈ R^MIt is identical.

The psnr of three kinds of computational methods of result of calculation second case is 118.7734, three kinds of sides of third type The psnr of method is 144.3529, MSE=2.38611e^-010.As can be seen that algorithm GR_OMP and QR_OMP tool given herein There is obviously advantage, especially has the characteristics that large-scale problem fast and effective.Such as the third situation above, it gives herein / 10th of time less than OMP algorithms of calculating of the algorithm gone out, such as table 2 (2-D Lena images (256 × 256)) and 3 (2- of table The case where D (1024 × 1024)) shown in：

Table 2

	IGR_OMP	OMP	IRLS	SP	GBP	IHT	CoSaMP
								CPUtime(s)	0.18	0.31	8.37	1.40	4.18	0.23	2.30
Psnr	22.06	22.08	23.63	21.04	22.77	15.66	19.77

Table 3

	IQR_OMP	OMP	IRLS	SP	GBP	IHT	CoSaMP
								CPUtime(s)	11.30	31.86	5809	243	2325	14.03	533
Psnr	45.90	45.90	48.70	46.34	47.83	46.79	45.27

(1) in table 2 256 × 256 Lena images, calculation matrix is random matrix, and sparse transformation dct transform is by column It is calculated.Wherein it is respectively：

IGR_OMP：It is increment OMP algorithm of the foundation given herein on the basis of Greville recursion；

IQR_OMP：It is increment OMP algorithm of the foundation given herein on QR decomposition bases；

OMP：It is OMP method of the foundation of classics in normal equation；

IRLS：Iteration weighted least-squares method (Iteratively Reweighted Least Square)；

SP：Subspace method for tracing；

GBP：Greedy base method for tracing (Greedy Basis Pursuit)；

IHT：Iterative hard thresholding method；

CoSaMP：Compression sampling match tracing method.

IQR_OMP, IGR_OMP algorithm are the method for the present invention in these methods, Hong Kong section that other methods are all The software that skill university Chengfu doctors Huo provide.

It is apparent that result of calculation shows that delta algorithm IGR_OMP algorithms presented and discussed herein and IQR_OMP algorithms have Superiority.

The present invention provides a kind of increment orthogonal matching pursuit methods of higher-dimension sparse vector reconstruct, implement the technology There are many method and approach of scheme, the above is only a preferred embodiment of the present invention, it is noted that for the art Those of ordinary skill for, various improvements and modifications may be made without departing from the principle of the present invention, these change Protection scope of the present invention is also should be regarded as into retouching.The available prior art of each component part being not known in the present embodiment adds To realize.

Claims (3)

1. a kind of higher-dimension sparse vector reconstructing method based on IQR_OMP, which is characterized in that include the following steps：

Step B1, input observation data vector y ∈ R^mWith dictionary matrix A=[a₁, a₂..., a_n]∈R^m×n, initial support collection is enabled to refer to MarkEnable initial residual vector r₀Equal to observation data vector y, i.e. r₀=y enables k=1, wherein R^mIndicate that m × 1 ties up real number Column vector, R^m×nIndicate that m × n ties up real number matrix, a_iIt represents i-th of m × 1 in dictionary matrix A and ties up real number column vector, i=1, 2…,n；

Step B2, ask in dictionary matrix A with -1 iteration residual vector r of kth_k-1Kth time iteration strongest correlation row j_k：

j_k∈argmax_j|(r_k-1, a_j) |, Ω_k=Ω_k-1∪{j_k}；

Wherein Ω_kFor kth time iteration supported collection index, Ω_k-1For -1 iteration supported collection index of kth, a_jFor in dictionary matrix A Jth row real number column vector；

Step B3, by the Q of -1 iteration of kth_k-1And R_k-1The Q of recurrence calculation kth time iteration_kAnd R_k；

Step B4, update kth time iteration residual vectorWherein q_kIt is kth time Iterative Matrix Q_k K column vector；

Step B5, judges whether to meet stopping criterion, i.e. the Euclid norm of kth time residual vector meets ‖ r_k‖₂≤ ε, if Meet, executes step B6；Otherwise k is updated to k+1, and repeats step B2 to step B4；

Step B6 solves minimization problemSolution x_k：

WhereinFor kth time iteration support index Ω_kCorresponding dictionary matrix,j₁, j₂..., j_k∈Ω_k,It is matrix R_kGeneralized inverse,It is matrix Q_kTransposition；

Step B7 exports sparse coefficient vector x：

Wherein x (i) indicates i-th of element in sparse coefficient vector x, x_k(i) vector x is indicated_kIn i-th of element.

2. according to the method described in claim 1, it is characterized in that, step B3 includes：If k=1, j_kFor j₁, thenQ₁The 1st of matrix is classified as Q₁(:, 1) or q₁It is equal toI.e. Step B4 is turned to,

If k ≠ 1, following calculate is executed successively：

q₁=a_jk-Q_k-1β_k,

q_k=q₁/α_k,

Q_k=[Q_k-1:q_k],It is wherein vectorialIt is the jth of dictionary matrix A_kA column vector, q_kIt is matrix Q K-th of column vector, β_kVector, α are replaced by what kth time executed temporarily_kNumerical value is replaced by what kth time executed temporarily.

3. according to the method described in claim 2, it is characterized in that, in step B6,Representation vector Euclid norm.