Decrease in Computational Load and Increase in Accuracy for Filtering of Random Signals

1. Introduction

1.1. Preliminaries

Although there is an extensive literature on studying methods of filtering of random signals, the problem we consider seems to be different from the known ones. The objective of this paper is to formulate and solve this problem under scenarios which relate to a reduction of computation complexity, for the case of large covariance matrices, and an increase in the filtering accuracy.

Let $\Omega$ be a set of outcomes in probability space $(\Omega ,\Sigma ,\mu )$ for which $\Sigma$ is a $\sigma$–algebra of measurable subsets of $\Omega$ and $\mu :\Sigma \to [0,1]$ is an associated probability measure. Let $\mathbf{x}\in L^2(\Omega ,{\mathbb{R}}^m)$ be a signal of interest and $\mathbf{y}\in L^2(\Omega ,{\mathbb{R}}^n)$ be an observable signal.. Here $L^2(\Omega ,{\mathbb{R}}^m)$ is the space of square-integrable functions defined on $\Omega$ with values in ${\mathbb{R}}^m$, i.e., such that ${\int }_{\Omega }{\sum }_{i=1}^m{\left[{\mathbf{x}}_{\left(i\right)}\left(\omega \right)\right]}^{2}d\mu \left(\omega \right)<\infty .$

Further, let us write $\mathbf{x}={[{\mathbf{x}}_{\left(1\right)},\dots ,{\mathbf{x}}_{\left(m\right)}]}^T$ where ${\mathbf{x}}_{\left(i\right)}\in L^2(\Omega ,\mathbb{R})$, for $i=1,\dots ,m$. We assume without loss of generality that all random vectors are of the zero mean. Let us denote

$$\begin{array}{c}\hfill {\parallel \mathbf{x}\parallel }_{\Omega }^2={\int }_{\Omega }\sum _{i=1}^m{\left[{\mathbf{x}}_{\left(i\right)}\left(\omega \right)\right]}^{2}d\mu \left(\omega \right),E_{xy}={\int }_{\Omega }\mathbf{x}\left(\omega \right){\left[\mathbf{y}\left(\omega \right)\right]}^{T}d\mu \left(\omega \right),\end{array}$$

(1)

where $E_{xy}$ is the covariance matrix formed from $\mathbf{x}$ and $\mathbf{y}$.

Let filter F be represented by $F\left(\mathbf{y}\right)=M\mathbf{y}$ where $M\in {\mathbb{R}}^{m\times n}$. Note that each matrix $M\in {\mathbb{R}}^{m\times n}$ defines a bounded linear transformation $\mathcal{M}:L^2(\Omega ,{\mathbb{R}}^n)\to L^2(\Omega ,{\mathbb{R}}^m)$. It is customary to write M rather then $\mathcal{M}$ since $\left[\mathcal{M}\right(\mathbf{y}\left)\right]\left(\omega \right)=M\left[\mathbf{y}\right(\omega \left)\right]$, for each $\omega \in \Omega$.

The known problem of finding optimal filter F is formulated as follows. Given large covariance matrices $E_{xy}$ and $E_{yy}$, find M that solves

$$\underset{M}{min}{∥\mathbf{x}-M\mathbf{y}∥}_{\Omega }^2.$$

(2)

The minimal Frobenius norm solution to (2) is given by (see, for example, [1,2,3,4,5])

$$M=E_{xy}{E}_{yy}^{†},$$

(3)

where $E_{yy}^{†}$ is the Moor-Penrose pseudo-inverse of $E_{yy}$.

1.2. Motivations

In 3, M contains $m\times n$ entries. In other words, the filter defined by 3 contains $m\times n$ parameters to optimize the associated accuracy.

Motivation 1: Increase in accuracy. Since M in 3 is optimal there is no other linear filter that can provide a better associated accuracy. At the same time, it may happen that the accuracy is not satisfactory. Then a logical way to increase the accuracy is an increase in the number of parameters to optimize. Since m is fixed then only n can be varied; n can be increased to, say, q where $q>n$, if $\mathbf{y}$ is replaced with another vector $\mathbf{v}\in L^2(\Omega ,{\mathbb{R}}^q)$ which is unknown. In this case, we need to determine both M and $\mathbf{v}$. In particular, $\mathbf{v}$ can be represented as $\mathbf{v}={[{\mathbf{v}}_1^T,\dots ,{\mathbf{v}}_p^T]}^T$ where, for $j=1,\dots ,p,$${\mathbf{v}}_j\in L^2(\Omega ,{\mathbb{R}}^{q_j})$ and $q=q_1+\dots +q_p$. Then the associated filter $F_p$ is represented by the equation

$$F_p({\mathbf{v}}_1,\dots {\mathbf{v}}_p)=M\mathbf{v}=[M_1,\dots ,M_p]{[{\mathbf{v}}_1^T,\dots ,{\mathbf{v}}_p^T]}^T=\sum _{j=1}^p{M}_j{\mathbf{v}}_j,$$

(4)

where, for $j=1,\dots ,p,$$M_j\in {\mathbb{R}}^{m\times q_j}$. We call p the degree of filter $F_p$.

Therefore, the preliminary statement of the problem considered here and in [6] is as follows. Given $E_{xv}$ and $E_{vv}$, find $M_1,\dots ,$$M_p$ and ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$ that solve

$$\underset{\begin{array}{c}{M}_1,\dots ,M_p\\ {\mathbf{v}}_1,\dots ,{\mathbf{v}}_p\end{array}}{min}{∥\mathbf{x}-\sum _{j=1}^p{M}_j{\mathbf{v}}_j∥}_{\Omega }^2.$$

(5)

Note, since $E_{xv}$ and $E_{vv}$ are given then, for $i,j=1,\dots ,p$, $E_{x{v}_j}$ and $E_{v_i{v}_j}$ are known as blocks of $E_{xv}$ and $E_{vv}$.

Motivation 2: Decrease in computational load. In 2 and 3, dimensions of $\mathbf{x}$ and $\mathbf{y}$, m and n, respectively, can be very large. In a number of important applied problems, such as those in biology, ecology, finance, sociology and medicine (see e.g., [7,8]), $m,n=O\left({10}^3\right)$ and greater. For example, measurements of gene expression can contain tens of hundreds of samples. Each sample can contain tens of thousands of genes. As a result, in such cases, associated covariance matrices matrices in 3 are very large. Further, in (4) and (5), dimension $q_j$ of ${\mathbf{v}}_j$ might be of the same order as that of observable signal $\mathbf{y}$. For instance, one of ${\mathbf{v}}_j$ can be equal to $\mathbf{y}$. Based on 2–(5), it is natural to predict that a solution of problem (5) contains $q\times q$ pseudo-inverse $E_{vv}^{†}$ and $m\times q$ pseudo-inverse $M^{†}$. It is indeed so, as shown in Section 4.1.6 and [6]. Matrices $E_{vv}^{†}$ and $M^{†}$ are usually larger than $E_{yy}^{†}$. Then computation of a large pseudo-inverse matrix leads to a quite slow computing procedure up to the cases when the computer runs out of computer memory. In particular, the associated well-known phenomenon of the “curse of dimensionality” [9] states that many problems become exponentially difficult in high dimensions.

Of course, in the case of high dimensional signals, some well known dimensionality reduction techniques can be used (as those, for example, in [2,3,5]) as an intermediate step for the filtering. At the same time, those techniques also involve large pseudo-inverse covariance matrices of the same size as that in 3. Additionally, in applying the dimensionality reduction procedure, intrinsic associated errors appear. The errors may involve, in particular, an undesirable distortion of the data. Thus, an application of the dimensionality reduction techniques to the problem under consideration does not allow us to avoid a computation of large pseudo-inverse covariance matrices. Therefore, we wish to develop techniques that reduce computation of large pseudo-inverses $E_{vv}^{†}$ and $M^{†}$ without involving any additional associated error.

The following example illustrates Motivation 2. Simulations in this example and also in 2 of Section 4.1.7 were run on a Dell Latitude 7400 Business Laptop (manufactured in 2020).

Example 1.

Let $E_{vv}^{†}\in {\mathbb{R}}^{q\times q}$ be a matrix whose entries are normally distributed with mean 0 and variance 1. In Table 1, we provide the time (in seconds) needed forMatlabto compute matrix $E_{vv}^{†}$ for different values of q.

We observe that computation of, e.g., $5000\times 5000$ pseudo-inverse matrix requires 63 seconds while computation of the five pseudo-inverse matrices of sizes $1000\times 1000$ requires about $2=5\times 0.4$ seconds. Thus, a procedure that would allow us to reduce the computation of $E_{vv}^{†}\in {\mathbb{R}}^{5000\times 5000}$ to a computation of, say, five pseudo-inverse matrices of sizes $1000\times 1000$ would be faster. At the same time, such a procedure would additionally contain associated matrix multiplications which are time consuming as well. This observation is further elaborated in Section 4.1.2 and 2 of Section 4.1.7.

2. Statement of the Problem

The solution of the problem in (5) is divided into two parts. In this paper, we consider the first part of the solution. It mainly concerns a technique for the decrease in a computational load associated with computation of $E_{vv}^{†}$. The second part of the solution of problem (5) is considered in [6]. It represents an extension of the technique considered here to the solution of the original problem in (5). More precisely, here, we consider

$$\underset{M_1,\dots ,M_p}{min}{∥\mathbf{x}-\sum _{j=1}^p{M}_j{\mathbf{v}}_j∥}_{\Omega }^2.$$

(6)

In terms of notation (4), the known minimal Frobenius norm solution of (6) is given by (see, e.g., [1,2,3,4,5])

$$M=M=E_{xv}{E}_{vv}^{†}.$$

(7)

As mentioned in Section 1.2, in (7) (and in (8) below), $q\times q$ matrix $E_{vv}$ might be very large. Then the pseudo-inverse $E_{vv}^{†}$ may be difficult to compute. In this case, formula in (7) is not useful in practice.

Therefore, the problem is as follows. Given large $E_{xv}$ and $E_{vv}$, find a computational procedure for solving (6) which significantly decreases the computational complexity for the $E_{vv}^{†}$ evaluation compared to that needed by the filter represented by (7).

While problem (6) arises in the solution provided in [6] it is also important in its own right. In particular, if $\mathbf{x}={[{\mathbf{x}}_1^T,\dots ,{\mathbf{x}}_p^T]}^T$ where, for $j=1,\dots ,p$, ${\mathbf{x}}_j\in L^2(\Omega ,{\mathbb{R}}^{m_j})$ and $m=m_1+\cdots +m_p$, then (6) represents the problem of finding the optimal multi-input multi-output (MIMO) filter. In this case, ${\mathbf{v}}_j$ and ${\mathbf{x}}_j$, for $j=1,\dots ,p$, are an input signal and output signal, respectively.

The problem in (6) can also be interpreted as the system identification problem [10,11,12,13]. In this case, ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$ are the inputs of the system, $M=[M_1,\dots ,M_p]$ is the input-output map of the system, and $\mathbf{x}$ is the system output. For example, in environmental monitoring, ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$ may represent random data coming from nodes measuring temperature, light or pressure variations, and $\mathbf{x}$ may represent a random signal obtained after a merging of the received data in order to estimate required data parameters within a prescribed accuracy. Similar scenarios occur, e.g., in target localization and tracking [14].

As mentioned above, in all problems considered in this paper, covariance matrices are assumed to be known. This assumption is a well-known and significant limitation in problems dealing with the optimal filtering of random signals. See [15,16,17,18] in this regard. The covariances can be estimated in various ways and particular techniques are given in many papers. We cite [19,20,21,22,23,24] as the examples. Estimation of covariance matrices is an important problem which is not considered in this paper.

3. Contribution and Novelty

The main conceptual novelty of the approach developed in this paper is in the technique that allows us to decrease the computational load needed for the numerical realization of the proposed filter. The procedure is based on the reduction of a large $E_{vv}$ matrix to a collection of smaller matrices. It is done so that the filter equation with large matrices is equivalently represented by the set of equations with smaller matrices. As a result, as shown in Section 4.1.7, the proposed p-th degree filter requires computation of p pseudo-inverse matrices with sizes smaller that the size of the large matrix $E_{vv}$. In this regard, see Section 4.1.2, and 1 and 2, for more details. The associated Matlab code is given in Section 4.1.7.

In 2, 3 and 4, we also show that the accuracy of the proposed filter increases with increases in the degree of the filter. Details are provided in Section 4.1 and more specifically, in Section 4.1.2. In particular, the algorithm for the optimal determination of $M_1,\dots ,$$M_p$ is given in Section 4.1.7. Note that the algorithm is easy to implement numerically. The error analysis associated with the filter under consideration is provided in Section 4.1.3, Section 4.1.5 and Section 4.1.8.

Remark 1.

Note that for every outcome $\omega \in \Omega$, realizations ${\mathbf{v}}_1\left(\omega \right),$$\dots ,{\mathbf{v}}_p\left(\omega \right)$ and $\mathbf{x}\left(\omega \right)$ of signals ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$ and $\mathbf{x}$, respectively, occur with certain probabilities. Thus, random signals ${\mathbf{v}}_j$, for $j=1,\dots ,p,$ and $\mathbf{x}$ are associated with infinite sets of realizations $V_j=\left\{{\mathbf{v}}_j\left(\omega \right)\right|\omega \in \Omega \}\subseteq {\mathbb{R}}^n$ and $X=\left\{\mathbf{x}\left(\omega \right)\right|\omega \in \Omega \}\subseteq {\mathbb{R}}^m$, respectively. For different ω, the values ${\mathbf{v}}_j\left(\omega \right)$, for $j=1,\dots ,p,$ and $\mathbf{x}\left(\omega \right)$ are, in general, different, and in many cases will span the entire spaces ${\mathbb{R}}_j^q$ and ${\mathbb{R}}^m$, respectively. Therefore, for each $\mathbf{x}$ and ${\mathbf{v}}_j$, for $j=1,\dots ,p$, filter given by (4) can be interpreted as map ${\mathbb{R}}^{q_1}\times \dots {\mathbb{R}}^{q_p}\to {\mathbb{R}}^m$ where probabilities are assigned to each column of ${\mathbb{R}}^{q_j}$ and ${\mathbb{R}}^m$. Importantly, the filter is invariant with respect to outcome $\omega \in \Omega$, and therefore, is the same for all ${\mathbf{v}}_1\left(\omega \right)\dots ,{\mathbf{v}}_p\left(\omega \right)$ and $\mathbf{x}\left(\omega \right)$ with different outcomes $\omega \in \Omega$.

4. Solution of the Problem

We wish that in (5) and (6), ${\sum }_{j=1}^p{M}_j{\mathbf{v}}_j$ would never been the zero vector. To this end, we need to extend the known definition of the linear independence of vectors as follows.

Let $M^{\left(1\right)}\in {\mathbb{R}}^{m\times q_1},...,M^{\left(p\right)}\in {\mathbb{R}}^{m\times q_p}$ be some matrices. For $j=1,\dots ,p$, let $\mathcal{N}\left(M^{\left(j\right)}\right)$ be a null space of matrix $M^{\left(j\right)}$.

A linear combination in the generalized sense of vectors ${\mathbf{v}}_1,\dots ,$${\mathbf{v}}_p$ is a vector

$$\widehat{\mathbf{v}}=M^{\left(1\right)}{\mathbf{v}}_1\left(\omega \right)+\dots +M^{\left(p\right)}{\mathbf{v}}_p\left(\omega \right).$$

We say that the linear combination in generalized sense is non-trivial if ${\mathbf{v}}_j\left(\omega \right)\notin \mathcal{N}\left(M^{\left(j\right)}\right)$, for at least one of $j=1,\dots ,p.$ Note that if $M^{\left(j\right)}\ne \mathbb{O}$ where $\mathbb{O}$ is the zero matrix, then it is still, of course, possible that ${\mathbf{v}}_j\left(\omega \right)\in \mathcal{N}\left(M^{\left(j\right)}\right)$. Therefore, condition ${\mathbf{v}}_j\left(\omega \right)\notin \mathcal{N}\left(M^{\left(j\right)}\right)$ is more general than condition $M^{\left(j\right)}\ne \mathbb{O}$.

Definition 1.

Random vectors ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$ are called linearly independent in the generalized sense if there is no non-trivial linear combination in the generalized sense of these vectors equal to the zero vector; in other words, if

$$M^{\left(1\right)}{\mathbf{v}}_1\left(\omega \right)+\dots +M^{\left(p\right)}{\mathbf{v}}_p\left(\omega \right)=\mathbb{O},$$

for almost all $\omega \in \Omega$, implies that ${\mathbf{v}}_j\left(\omega \right)\in \mathcal{N}\left(M^{\left(j\right)}\right)$, for each $j=1,\dots ,p$, and almost all $\omega \in \Omega$.

All random vectors considered below are assumed to be linearly independent in the generalized sense.

4.1. Solution of Problem in (6)

In (7), matrix M is the minimum Frobenius norm solution of the equation

$$M{E}_{vv}=E_{xv}.$$

(8)

Here, we exploit the specific structure of the original system (8) to develop a special block elimination procedure that allows us to reduce the original system of equations (8) to a collection of independent smaller subsystems. Their solution implies less computational load than that needed for (7).

To this end, let us represent $E_{vv}$ and $E_{xv}$ in blocked forms as follows:

$$E_{vv}=\left[\begin{array}{ccc}{E}_{11}& \dots & E_{1p}\\ ⋮& ⋮& ⋮\\ E_{p1}& \dots & E_{pp}\end{array}\right]\mathrm{and}{E}_{xv}=[E_1\dots E_p],$$

(9)

where $E_{ij}=E_{v_i{v}_j}\in {\mathbb{R}}^{q_j\times q_j}$ and $E_j=E_{x{v}_j}\in {\mathbb{R}}^{m\times q_j}$, for $j,i=1,\dots ,p$. Then (8) implies

$$[M_1\dots M_p]\left[\begin{array}{ccc}{E}_{11}& \dots & E_{1p}\\ ⋮& ⋮& ⋮\\ E_{p1}& \dots & E_{pp}\end{array}\right]=[E_1\dots E_p].$$

(10)

Thus, the solution of the problem in (8) is equivalent to the solution of equation (10).

We provide a specific solution of equation (10) in the following way. First, in Section 4.1.1, a generic basis for the solution is given. In Section 4.1.4, the solution of equation (10) is considered for $p=3$. This is a preliminary step for the solution of equation (10) for an arbitrary p which is provided in Section 4.1.6.

4.1.1. Generic Basis for the Solution of Problem (8): Case $p=2$ in (10). Second Degree Filter

To provide a generic basis for the solution of equation (10), let us consider the case $p=2$ in (4) and (6), i.e., the second degree filter. Then (10) becomes

$$\left[M_1{M}_2\right]\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& E_{22}\end{array}\right]=\left[E_1{E}_2\right].$$

(11)

We denote

$$E_{22}^{\left(1\right)}=E_{22}-E_{21}{E}_{11}^{†}{E}_{12}\mathrm{and}{E}_2^{\left(1\right)}=E_2-E_1{E}_{11}^{†}{E}_{12}.$$

(12)

Recall, it has been shown in [4] that for any random vectors g and h,

$$E_{gh}=E_{gh}{E}_{hh}^{†}{E}_{hh}.$$

(13)

Lemma 1.

Let $M_1$ and $M_2$ satisfy equation (11). Then (11) decomposes into equations

$$M_1{E}_{11}=E_1-M_2{E}_{21}\mathrm{and}{M}_2{E}_{22}^{\left(1\right)}=E_2^{\left(1\right)}$$

(14)

Proof. 

Let us denote ${\overline{M}}_1=M_1+M_2{E}_{21}{E}_{11}^{†}$. Then

$$\left[{\overline{M}}_1{M}_2\right]=\left[M_1{M}_2\right]\left[\begin{array}{cc}I& \mathbb{O}\\ E_{21}{E}_{11}^{†}& I\end{array}\right]$$

and

$$\left[M_1{M}_2\right]=\left[{\overline{M}}_1{M}_2\right]\left[\begin{array}{cc}I& \mathbb{O}\\ -E_{21}{E}_{11}^{†}& I\end{array}\right],$$

(15)

where by (13),

$$E_{21}{E}_{11}^{†}{E}_{11}=E_{21}.$$

(16)

Thus, (11) is represented as

$$\left[{\overline{M}}_1{M}_2\right]\left[\begin{array}{cc}I& \mathbb{O}\\ -E_{21}{E}_{11}^{†}& I\end{array}\right]\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& E_{22}\end{array}\right]=\left[E_1{E}_2\right]$$

which implies

$$\left[{\overline{M}}_1{M}_2\right]\left[\begin{array}{cc}I& \mathbb{O}\\ -E_{21}{E}_{11}^{†}& I\end{array}\right]\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& E_{22}\end{array}\right]\left[\begin{array}{cc}I& -E_{11}^{†}{E}_{12}\\ \mathbb{O}& I\end{array}\right]$$

$$=\left[E_1{E}_2\right]\left[\begin{array}{cc}I& -E_{11}^{†}{E}_{12}\\ \mathbb{O}& I\end{array}\right].$$

(17)

On the basis of (17), we see that (8), for $p=2$, reduces to the equation

$$\left[{\overline{M}}_1{M}_2\right]\left[\begin{array}{cc}{E}_{11}& \mathbb{O}\\ \mathbb{O}& E_{22}^{\left(1\right)}\end{array}\right]=\left[E_1{E}_2^{\left(1\right)}\right].$$

(18)

The latter implies 14. ▪

Remark 2.

The above proof is based on equation (15) which implies the equivalence of equations (11) and 14. In turn, (1) allows us to equivalently obtain 14 as follows. Let us multiply (11) from the right by $N_{21}=\left[\begin{array}{cc}I& -{E_{11}}^{†}{E}_{12}\\ \mathbb{O}& I\end{array}\right]$. Then

$$\left[M_1{M}_2\right]\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& E_{22}\end{array}\right]N_{21}=\left[E_1{E}_2\right]N_{21},$$

(19)

i.e.,

$$\left[M_1{M}_2\right]\left[\begin{array}{cc}{E}_{11}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}\end{array}\right]=\left[E_1{E}_2^{\left(1\right)}\right].$$

(20)

Then 14 follows from (20).

Theorem 1.

For $p=2,$ the minimal Frobenius norm solution of problem (8) is represented by

$$M_2={\tilde{M}}_2=E_2^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†},M_1={\tilde{M}}_1=(E_1-{\tilde{M}}_2{E}_{21})E_{11}^{†}$$

(21)

Proof. 

By (1), matrices $M_1$ and $M_2$ that solve (11) satisfy the equations in 14. By [2,3,17,18], their minimal Frobenius norm solutions are given by (21). ▪

It follows from (21) that the second degree filter requires computation of two pseudo-inverse matrices, $E_{11}^{†}\in {\mathbb{R}}^{q_1\times q_1}$ and ${\left(E_{22}^{\left(1\right)}\right)}^{†}\in {\mathbb{R}}^{q_2\times q_2}$, of sizes that are smaller than the size of $E_{vv}^{†}\in {\mathbb{R}}^{q\times q}$ (if, of course, $q_1\ne 0$ and $q_2\ne 0$).

4.1.2. Decrease in Computational Load

Here, we wish to compare the computational load needed for a numerical realization of the second degree filter given, for $p=2$, by (4) and (21), and that of the filter represented by (7).

By [25] (p. 254), the product of $m\times n$ and $n\times q$ matrices requires, for large $m,n$ and q, $2mnq$ flops. For large q, the Golub-Reinsch SVD used for computation of the $q\times q$ matrix pseudo-inverse requires $21q^3$ flops (see [25], p. 254).

In (12), for $q_1=q_2=q/2$, the evaluation of $E_{22}^{\left(1\right)}$ implies the evaluation of $E_{11}^{†}$, matrix multiplications in $E_{21}{E}_{11}^{†}{E}_{12}$ and matrix subtraction. These operations require $21{(q/2)}^3$ flops, $2{(q/2)}^3+2{(q/2)}^3$ flops and $q^2/4$ flops, respectively. Similarly, computation of $E_2^{\left(1\right)}$ requires $2m{q}^2/4+2q^3/8+mq/2$ flops. In (21), computation of $M_2$ and $M_1$ imply $21q^3/8+2m{q}^2/4$ flops and $2m{q}^2/4+2m{q}^2/4+mq/2$ flops, respectively.

Let us denote by $C\left(A\right)$ a number of flops required to compute an operation A. In total, the second degree filter requires $C\left(E_{22}^{\left(1\right)}\right)+C\left(E_2^{\left(1\right)}\right)+C\left(M_2\right)+C\left(M_1\right)$ flops where

$$C\left(E_{22}^{\left(1\right)}\right)=q^2/4+25q^3/8,$$

$$C\left(E_2^{\left(1\right)}\right)=mq/2+2m{q}^2/4+2q^3/8,$$

$$C\left(M_2\right)=2m{q}^2/4+21q^3/8,$$

$$C\left(M_1\right)=mq/2+m{q}^2.$$

Thus, the number of flops needed for the numerical realization of the second degree filter is equal to $mq+q^2/4+2m{q}^2+6q^3$.

The computational load required by the filter given by (7) is $2m{q}^2+21q^3$ flops. It is customary to choose $m\le q$. In this case, clearly, $mq+q^2/4+2m{q}^2+6q^3<2m{q}^2+21q^3$. In particular, if $m=q$ then the latter inequality becomes $5q^2/4+8q^3<23q^3.$ That is, for sufficiently large m and q, the proposed second degree filter requires, roughly, up to $2.8$ times less flops than that by the filter given by (7). In simulations by Matlab provided in 2 below, the second degree filter takes about two times less seconds than the filter given by (7). A command that would allow us to compute a number of flops is not available in Matlab.

In the following Section 4.1.4 and Section 4.1.6, the procedure described in Section 4.1.1 is extended to the case of the filter of higher degrees. For given q, this will imply, obviously, smaller sizes of blocks of matrices $E_{vv}$ and $E_{xv}$ than those for the second degree filter, and as a result, a further decrease in the associated computational load.

4.1.3. Error Associated with the Second Degree Filter Determined by (4)

Now, we wish to show that the error associated with the second degree filter determined, for $p=2$, by (4) and (21) is less than the error associated with the filter determined by 3.

Let us denote

$$ϵ=\underset{M}{min}{∥\mathbf{x}-M\mathbf{y}∥}_{\Omega }^2,{\epsilon }_2=\underset{M_1,M_2}{min}{∥\mathbf{x}-\left[M_1{M}_2\right]\left[\begin{array}{c}{\mathbf{v}}_1\\ {\mathbf{v}}_2\end{array}\right]∥}_{\Omega }^2.$$

(22)

The Frobenius norm of matrix M is denoted by $\parallel M\parallel$.

Theorem 2.

Let ${\mathbf{v}}_1$ and ${\mathbf{v}}_2$ be such that

$$\sum _{j=1}^2{∥E_j^{(j-1)}{\left[{\left(E_{jj}^{(j-1)}\right)}^{1/2}\right]}^{†}∥}^2>{\parallel E_{xy}{E}_{yy}^{1/2}\parallel }^2,$$

(23)

where $E_1=E_1^{\left(0\right)}$ and $E_{11}=E_{11}^{\left(0\right)}$. Then

$${\epsilon }_2<ϵ.$$

(24)

Proof. 

It is known (see, for example, [1,2,3,4,5]) that, for M determined by 3, and for ${\tilde{M}}_1$ and ${\tilde{M}}_2$ determined by (21),

$$ϵ=\parallel E_{xx}^{1/2}{\parallel }^2-\mathrm{tr}\left\{E_{xy}{E}_{yy}^{†}{E}_{yx}\right\},{\epsilon }_2={\parallel E_{xx}^{1/2}\parallel }^2-\mathrm{tr}\left\{E_{xv}{E}_{vv}^{†}{E}_{vx}\right\},$$

(25)

where $\mathbf{v}=\left[\begin{array}{c}{\mathbf{v}}_1\\ {\mathbf{v}}_2\end{array}\right]$ and, bearing in mind (2), $\left[{\tilde{M}}_1{\tilde{M}}_2\right]=E_{xv}{E}_{vv}^{†}$. Therefore, by (21),

$$\mathrm{tr}\left\{E_{xv}{E}_{vv}^{†}{E}_{vx}\right\}=\mathrm{tr}\left\{\left[{\tilde{M}}_1{\tilde{M}}_2\right]E_{vx}\right\}=\mathrm{tr}\{{\tilde{M}}_1{E}_1^T+{\tilde{M}}_2{E}_2^T\}$$

$$=\mathrm{tr}\{(E_1{E}_{11}^{†}-{\tilde{M}}_2{E}_{21}{E}_{11}^{†})E_1^T+{\tilde{M}}_2{E}_2^T\}$$

$$=\mathrm{tr}\left\{E_1{E}_{11}^{†}{E}_1^T\right\}+\mathrm{tr}\left\{{\tilde{M}}_2(E_2^T-E_{21}{E}_{11}^{†}{E}_1^T)\right\}$$

$$=\mathrm{tr}\left\{E_1{E}_{11}^{†}{E}_1^T\right\}$$

$$+\mathrm{tr}\left\{(E_2-E_1{E}_{11}^{†}{E}_{12}){(E_{22}-E_{21}{E}_{11}^{†}{E}_{12})}^{†}(E_2^T-E_{21}{E}_{11}^{†}{E}_1^T)\right\}$$

$$=\parallel E_1{E}_{11}^{1/2}{\parallel }^2+{\parallel E_2^{\left(1\right)}{\left[{\left(E_{22}^{\left(1\right)}\right)}^{1/2}\right]}^{†}\parallel }^2.$$

(26)

Then

$${\epsilon }_2=\parallel E_{xx}^{1/2}{\parallel }^2-(\parallel E_1{E}_{11}^{1/2}{\parallel }^2+\parallel E_2^{\left(1\right)}{\left[{\left(E_{22}^{\left(1\right)}\right)}^{1/2}\right]}^{†}{\parallel }^2).$$

(27)

At the same time,

$$ϵ=\parallel E_{xx}^{1/2}{\parallel }^2-{\parallel E_{xy}{E}_{yy}^{1/2}\parallel }^2.$$

(28)

As a result, (24) follows from (23), (27) and (28). ▪

Remark 3.

The condition in (23) can practically be satisfied by, in particular, a suitable choice of ${\mathbf{v}}_1$ and ${\mathbf{v}}_2$, and/or their dimensions. The following 1 demonstrates such a particular choice of ${\mathbf{v}}_1$.

Corollary 1.

Let ${\mathbf{v}}_1=\mathbf{y}$ and ${\mathbf{v}}_2$ be a nonzero signal. Then (24) is always true.

Proof. 

The proof follows directly from (23), (27) and (28). ▪

In general, the proposed solution of equation ((10)) is based on a combination of (1) and (2) with the extension of the idea of Gaussian elimination [25] to the case of matrices with specific blocks $E_{ij}$ and $E_j$. This allows us to build the special procedure considered below in Section 4.1.4 and Section 4.1.6.

4.1.4. Solution of Equation (10) for $p=3$

To clarify the solution device of equation (10) for an arbitrary p, let us now consider a particular case with $p=3$. The solution is based on the development of the procedure represented by (19) and (20).

For $p=3$, equation (10) is as follows:

$$\left[M_1{M}_2{M}_3\right]\left[\begin{array}{ccc}{E}_{11}& E_{12}& E_{13}\\ E_{21}& E_{22}& E_{23}\\ E_{31}& E_{32}& E_{33}\end{array}\right]=\left[E_1{E}_2{E}_3\right].$$

(29)

First, consider a reduction of (29) to the equation with the block-lower triangular matrix.

Step 1. By (13), for $i=2,3$, $E_{i1}=E_{i1}{E}_{11}^{†}{E}_{11}$. Define

$$N_{3,1}=\left[\begin{array}{ccc}I& -{E_{11}}^{†}{E}_{12}& -{E_{11}}^{†}{E}_{13}\\ \mathbb{O}& I& \mathbb{O}\\ \mathbb{O}& \mathbb{O}& I\end{array}\right].$$

Now, in (29), for $i=1,2,3$ and $j=2,3$, update block-matrices $E_{ij}$ and $E_j$ to $E_{ij}^{\left(1\right)}$ and $E_j^{\left(1\right)}$, respectively, so that

$$\left[\begin{array}{ccc}{E}_{11}& \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& E_{23}^{\left(1\right)}\\ E_{31}& E_{32}^{\left(1\right)}& E_{33}^{\left(1\right)}\end{array}\right]=\left[\begin{array}{ccc}{E}_{11}& E_{12}& E_{13}\\ E_{21}& E_{22}& E_{23}\\ E_{31}& E_{32}& E_{33}\end{array}\right]N_{3,1}$$

(30)

and

$$\left[E_1{E}_2^{\left(1\right)}{E}_3^{\left(1\right)}\right]=\left[E_1{E}_2{E}_3\right]N_{3,1}$$

(31)

where $E_{kj}^{\left(1\right)}=E_{kj}-E_{k1}{E}_{11}^{†}{E}_{1j}$, for $k=2,3$ and $j=2,3$, and $E_j^{\left(1\right)}=E_j-E_1{E}_{11}^{†}{E}_{1j}$, for $j=2,3$. Note that $E_{12}^{\left(1\right)}=E_{13}^{\left(1\right)}=\mathbb{O}$ by (13).

Step 2. Set

$$N_{3,2}=\left[\begin{array}{ccc}I& \mathbb{O}& \mathbb{O}\\ \mathbb{O}& I& -{{E_{22}}^{\left(1\right)}}^{†}{E}_{23}^{\left(1\right)}\\ \mathbb{O}& \mathbb{O}& I\end{array}\right].$$

By (13), $E_{32}^{\left(1\right)}$$=E_{32}^{\left(1\right)}{E_{22}^{\left(1\right)}}^{†}{E}_{22}^{\left(1\right)}$. Taking the latter into account, update block-matrices $E_{ij}^{\left(1\right)}$ and $E_j^{\left(1\right)}$ in the LHS of (31) and 31, for $i=2,3$ and $j=3$, so that

$$\left[\begin{array}{ccc}{E}_{11}& \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \mathbb{O}\\ E_{31}& E_{32}^{\left(1\right)}& E_{33}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{ccc}{E}_{11}& \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& E_{23}^{\left(1\right)}\\ E_{31}& E_{32}^{\left(1\right)}& E_{33}^{\left(1\right)}\end{array}\right]N_{3,2}$$

(32)

and

$$\left[E_1{E}_2^{\left(1\right)}{E}_3^{\left(2\right)}\right]=\left[E_1{E}_2^{\left(1\right)}{E}_3^{\left(1\right)}\right]N_{3,2}$$

(33)

where $E_{33}^{\left(2\right)}=E_{33}^{\left(1\right)}-E_{32}^{\left(1\right)}{E_{22}^{\left(1\right)}}^{†}{E}_{23}^{\left(1\right)}$ and $E_3^{\left(2\right)}=E_3^{\left(1\right)}-E_2^{\left(1\right)}{E_{22}^{\left(1\right)}}^{†}{E}_{23}^{\left(1\right)}$.

As a result, we finally obtain the updated matrix equation with the block-lower triangular matrix

$$\left[M_1{M}_2{M}_3\right]\left[\begin{array}{ccc}{E}_{11}& \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \mathbb{O}\\ E_{31}& E_{32}^{\left(1\right)}& E_{33}^{\left(2\right)}\end{array}\right]=\left[E_1{E}_2^{\left(1\right)}{E}_3^{\left(2\right)}\right].$$

(34)

Solution of equation (34). Using a back substitution, equation (34) is represented by the three equations as follows:

$$M_3{E}_{33}^{\left(2\right)}=E_3^{\left(2\right)},$$

(35)

$$M_2{E}_{22}^{\left(1\right)}+M_3{E}_{32}^{\left(1\right)}=E_2^{\left(1\right)},$$

(36)

$$M_1{E}_{11}+M_2{E}_{21}+M_3{E}_{31}=E_1,$$

(37)

where $E_{11},E_{21},E_{31}$ and $E_1$ are the original blocks in (29).

It follows from (35)–(37) that the solution of equation (10), for $p=3$, is represented by the minimal Frobenius norm solutions to the problems

$$\underset{M_3}{min}{\parallel E_3^{\left(2\right)}-M_3{E}_{33}^{\left(2\right)}\parallel }^2,$$

(38)

$$\underset{M_2}{min}{\parallel (E_2^{\left(1\right)}-{\tilde{M}}_3{E}_{32}^{\left(1\right)})-M_2{E}_{22}^{\left(1\right)}\parallel }^2,$$

(39)

$$\underset{M_1}{min}{\parallel (E_1-\sum _{j=2}^3{\tilde{M}}_j{E}_{j1})-M_1{E}_{11}\parallel }^2,$$

(40)

where ${\tilde{M}}_3$ and ${\tilde{M}}_2$ are solutions to the problems in (38) and (39), respectively. Matrices ${\tilde{M}}_j$, for $j=1,2,3$, that solve (38)-(40) are as follows:

$${\tilde{M}}_3=E_3^{\left(2\right)}{\left[E_{33}^{\left(2\right)}\right]}^{†},$$

(41)

$${\tilde{M}}_2=(E_2^{\left(1\right)}-{\tilde{M}}_3{E}_{32}^{\left(1\right)}){\left(E_{22}^{\left(1\right)}\right)}^{†},$$

(42)

$${\tilde{M}}_1=(E_1-\sum _{j=2}^3{\tilde{M}}_j{E}_{j1}){\left(E_{11}\right)}^{†}.$$

(43)

Thus, the third degree filter requires computation of three pseudo-inverse matrices, $E_{11}^{†}\in {\mathbb{R}}^{q_1\times q_1}$, ${\left(E_{22}^{\left(1\right)}\right)}^{†}\in {\mathbb{R}}^{q_2\times q_2}$ and ${\left[E_{33}^{\left(2\right)}\right]}^{†}\in {\mathbb{R}}^{q_3\times q_3}$.

4.1.5. Error Associated with the Third Degree Filter Determined by (41)–(43)

We wish to show that the error associated with the filter determined by (41)–(43), is less than that associated with the filter determined by 3, $ϵ$.

To this end, let us denote

$${\epsilon }_3=\underset{M_1,M_2,M_3}{min}{∥\mathbf{x}-\left[M_1{M}_2{M}_3\right]\mathbf{v}∥}_{\Omega }^2,$$

(44)

where $\mathbf{v}={[{\mathbf{v}}_1^T,{\mathbf{v}}_2^T,{\mathbf{v}}_3^T]}^T$. As before, $E_1=E_1^{\left(0\right)}$ and $E_{11}=E_{11}^{\left(0\right)}$.

Theorem 3.

Let ${\mathbf{v}}_1$, ${\mathbf{v}}_2$ and ${\mathbf{v}}_3$ be such that

$$\sum _{j=1}^3{∥E_j^{(j-1)}{\left[{\left(E_{jj}^{(j-1)}\right)}^{1/2}\right]}^{†}∥}^2>{\parallel E_{xy}{E}_{yy}^{1/2}\parallel }^2.$$

(45)

Then

$${\epsilon }_3<ϵ.$$

(46)

Proof. 

We have

$${\epsilon }_3={\parallel E_{xx}^{1/2}\parallel }^2-\mathrm{tr}\left\{E_{xv}{E}_{vv}^{†}{E}_{vx}\right\},$$

(47)

and, bearing in mind an obvious extension of Remark (2), for $p=3$, $\left[{\tilde{M}}_1{\tilde{M}}_2{\tilde{M}}_3\right]=E_{xv}{E}_{vv}^{†}$. Therefore, by ((41))–((43)),

$$\begin{array}{c}\mathrm{tr}\left\{E_{xv}{E}_{vv}^{†}{E}_{vx}\right\}\hfill \\ =\mathrm{tr}\left\{\left[{\tilde{M}}_1{\tilde{M}}_2{\tilde{M}}_3\right]E_{vx}\right\}\hfill \\ =\mathrm{tr}\{{\tilde{M}}_1{E}_1^T+{\tilde{M}}_2{E}_2^T+{\tilde{M}}_3{E}_3^T\}\hfill \\ =\mathrm{tr}\{(E_1{E}_{11}^{†}-[{\tilde{M}}_2{E}_{21}+{\tilde{M}}_3{E}_{31}]E_{11}^{†})E_1^T+{\tilde{M}}_2{E}_2^T+{\tilde{M}}_3{E}_3^T\}\hfill \\ =\mathrm{tr}\left\{E_1{E}_{11}^{†}{E}_1^T\right\}+\mathrm{tr}\{{\tilde{M}}_2(E_2^T-E_{21}{E}_{11}^{†}{E}_1^T)\hfill \\ +{\tilde{M}}_3(E_3^T-E_{31}{E}_{11}^{†}{E}_1^T)\}\hfill \\ =\mathrm{tr}\left\{E_1{E}_{11}^{†}{E}_1^T\right\}+\mathrm{tr}\{(E_2^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†}-{\tilde{M}}_3{E}_{32}^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†}){\left(E_2^{\left(1\right)}\right)}^T\hfill \\ +{\tilde{M}}_3(E_3^T-E_{31}{E}_{11}^{†}{E}_1^T)\}\hfill \\ =\mathrm{tr}\{E_1{E}_{11}^{†}{E}_1^T+E_2^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†}{\left(E_2^{\left(1\right)}\right)}^T\hfill \\ +{\tilde{M}}_3[E_3^T-E_{31}{E}_{11}^{†}{E}_1^T-E_{32}^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†}+{\left(E_2^{\left(1\right)}\right)}^T]\}\hfill \\ =\parallel E_1{E}_{11}^{1/2}{\parallel }^2+{\parallel E_2{E}_{22}^{1/2}\parallel }^2+\mathrm{tr}\left\{E_3^{\left(2\right)}{\left(E_{33}^{\left(2\right)}\right)}^{†}{\left(E_3^{\left(2\right)}\right)}^T\right\}\hfill \\ =\parallel E_1{E}_{11}^{1/2}{\parallel }^2+\parallel E_2^{\left(1\right)}{\left[{\left(E_{22}^{\left(1\right)}\right)}^{1/2}\right]}^{†}{\parallel }^2+{\parallel E_3^{\left(2\right)}{\left[{\left(E_{33}^{\left(2\right)}\right)}^{1/2}\right]}^{†}\parallel }^2.\hfill \end{array}$$

(48)

Thus, (45)–(46) follows from $ϵ$ in (22), (47) and (48). ▪

Corollary 2.

If ${\mathbf{v}}_1=\mathbf{y}$, and at least one of ${\mathbf{v}}_2$ or ${\mathbf{v}}_3$ is a nonzero vector then (46) is always true.

Proof. 

The proof follows from (45), (48) and (28). ▪

4.1.6. Solution of Equation (10) for Arbitrary p

Here, we extend the procedure considered in (Section 4.1.4) to the case of an arbitrary p. First, we reduce the equation in (10) to an equation with the block-lower triangular matrix. We do it by the following steps where each step is justified by an obvious extension of (1).

Step 1. First, in (10), for $k=1,\dots ,p$ and $j=2,\dots ,p$, we wish to update blocks $E_{kj}$ and $E_k$ to $E_{kj}^{\left(1\right)}$ and $E_k^{\left(1\right)}$, respectively. To this end, define

$$N_{p,1}=\left[\begin{array}{cccc}I& -{E_{11}}^{†}{E}_{12}& \dots & -{E_{11}}^{†}{E}_{1p}\\ \mathbb{O}& I& \dots & \mathbb{O}\\ \dots & \dots & \dots & \dots \\ \mathbb{O}& \dots & I& \mathbb{O}\\ \mathbb{O}& \dots & \mathbb{O}& I\end{array}\right].$$

Bearing in mind that by (13), $E_{i1}=E_{i1}{E}_{11}^{†}{E}_{11}$, for $i=2,\dots ,p$, we obtain

$$\left[\begin{array}{cccc}{E}_{11}& \mathbb{O}& \dots & \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \dots & E_{2p}^{\left(1\right)}\\ \dots & \dots & \dots & \dots \\ E_{p-1,1}& E_{p-1,2}^{\left(1\right)}& \dots & E_{p-1,p}^{\left(1\right)}\\ E_{p1}& E_{p2}^{\left(1\right)}& \dots & E_{pp}^{\left(1\right)}\end{array}\right]$$

$$=\left[\begin{array}{cccc}{E}_{11}& E_{12}& \dots & E_{1p}\\ E_{21}& E_{22}& \dots & E_{2p}\\ \dots & \dots & \dots & \dots \\ E_{p-1,1}& E_{p-1,2}& \dots & E_{p-1,p}\\ E_{p1}& E_{p2}& \dots & E_{pp}\end{array}\right]N_{p,1}$$

(49)

and

$$[E_1{E}_2^{\left(1\right)}\dots E_p^{\left(1\right)}]=[E_1{E}_2\dots E_p]N_{p,1},$$

(50)

where, for $k=1,\dots ,p$ and $j=2,\dots ,p$,

$$E_{kj}^{\left(1\right)}=E_{kj}-E_{k1}{E}_{11}^{†}{E}_{1j}$$

(51)

and

$$E_j^{\left(1\right)}=E_j-E_1{E}_{11}^{†}{E}_{1j}.$$

(52)

Step 2. Now we wish to update blocks $E_{kj}^{\left(1\right)}$ and $E_k^{\left(1\right)}$, for $k=2,\dots ,p$ and $j=3,\dots ,p$, of the matrices in the LHS of (49) and (50), respectively.

Set

$$N_{p,2}=\left[\begin{array}{ccccc}I& \mathbb{O}& \mathbb{O}& \dots & \mathbb{O}\\ \mathbb{O}& I& -{{E_{22}}^{\left(1\right)}}^{†}{E}_{23}^{\left(1\right)}& \dots & -{{E_{22}}^{\left(1\right)}}^{†}{E}_{2p}^{\left(1\right)}\\ \mathbb{O}& \mathbb{O}& I& \dots & \mathbb{O}\\ \dots & \dots & \dots & \dots & \dots \\ \mathbb{O}& \mathbb{O}& \mathbb{O}& \dots & I\end{array}\right],$$

Because by (13), $E_{i2}^{\left(1\right)}=E_{i2}^{\left(1\right)}{E_{22}^{\left(1\right)}}^{†}{E}_{22}^{\left(1\right)}$, for $i=3,\dots ,p$, we have

$$\left[\begin{array}{ccccc}{E}_{11}& \mathbb{O}& \mathbb{O}& \dots & \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \mathbb{O}& \dots & \mathbb{O}\\ E_{31}& E_{32}^{\left(1\right)}& E_{33}^{\left(2\right)}& \dots & E_{3p}^{\left(2\right)}\\ \dots & \dots & \dots & \dots & \dots \\ E_{p1}& E_{p2}^{\left(1\right)}& E_{p3}^{\left(2\right)}& \dots & E_{pp}^{\left(2\right)}\end{array}\right]$$

$$=\left[\begin{array}{ccccc}{E}_{11}& \mathbb{O}& \mathbb{O}& \dots & \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& E_{23}^{\left(1\right)}& \dots & E_{2p}^{\left(1\right)}\\ E_{31}& E_{32}^{\left(1\right)}& E_{33}^{\left(1\right)}& \dots & E_{3p}^{\left(1\right)}\\ \dots & \dots & \dots & \dots & \dots \\ E_{p1}& E_{p2}^{\left(1\right)}& E_{p3}^{\left(1\right)}& \dots & E_{pp}^{\left(1\right)}\end{array}\right]N_{p,2}$$

(53)

and

$$[E_1{E}_2^{\left(1\right)}{E}_3^{\left(2\right)}\dots E_p^{\left(2\right)}]=[E_1{E}_2^{\left(1\right)}{E}_3^{\left(1\right)}\dots E_p^{\left(1\right)}]N_{p,2},$$

(54)

where, for $k=2,\dots ,p$ and $j=3,\dots ,p$,

$$E_{kj}^{\left(2\right)}=E_{kj}^{\left(1\right)}-E_{k2}^{\left(1\right)}{E_{22}^{\left(1\right)}}^{†}{E}_{2j}^{\left(1\right)}$$

(55)

and

$$E_j^{\left(2\right)}=E_j^{\left(1\right)}-E_2^{\left(1\right)}{E_{22}^{\left(1\right)}}^{†}{E}_{2j}^{\left(1\right)}.$$

(56)

We continue this updating procedure up to the final Step $(p-1)$. In Step $(p-2)$, the following matrices are obtained:

$$\left[\begin{array}{ccccc}{E}_{11}& \mathbb{O}& \dots & \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \dots & \mathbb{O}& \mathbb{O}\\ \dots & \dots & \dots & \dots & \dots \\ E_{p-1,1}& E_{p-1,2}^{\left(1\right)}& \dots & E_{p-1,p-1}^{(p-2)}& E_{p-1,p}^{(p-2)}\\ E_{p1}& E_{p2}^{\left(1\right)}& \dots & E_{p,p-1}^{(p-2)}& E_{pp}^{(p-2)}\end{array}\right]$$

(57)

and

$$[E_1{E}_2^{\left(1\right)}\dots E_{p-2}^{(p-3)}{E}_{p-1}^{(p-2)}{E}_p^{(p-2)}].$$

(58)

Step $(p-1)$. By (13),

$$E_{p-1,p}^{(p-2)}=E_{p-1,p}^{(p-2)}{\left(E_{p-1,p-1}^{\left(2\right)}\right)}^{†}{E}_{p-1,p-1}^{(p-2)}.$$

In 57 and 58, we now update only two blocks $E_{p-1,p}^{(p-2)}$ and $E_{pp}^{(p-2)}$, and the single block $E_p^{(p-2)}$. To this end, define

$$N_{p,p-1}=\left[\begin{array}{ccccc}I& \mathbb{O}& \dots & \mathbb{O}& \mathbb{O}\\ \mathbb{O}& I& \dots & \mathbb{O}& \mathbb{O}\\ \dots & \dots & \dots & \dots & \dots \\ \mathbb{O}& \mathbb{O}& \dots & I& -{E_{p-1,p-1}^{(p-2)}}^{†}{E}_{p-1,p}^{(p-2)}\\ \mathbb{O}& \mathbb{O}& \dots & \mathbb{O}& I\end{array}\right],$$

and obtain

$$\left[\begin{array}{ccccc}{E}_{11}& \mathbb{O}& \dots & \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \dots & \mathbb{O}& \mathbb{O}\\ \dots & \dots & \dots & \dots & \dots \\ E_{p-1,1}& E_{p-1,2}^{\left(1\right)}& \dots & E_{p-1,p-1}^{(p-2)}& \mathbb{O}\\ E_{p1}& E_{p2}^{\left(1\right)}& \dots & E_{p,p-1}^{(p-2)}& E_{pp}^{(p-1)}\end{array}\right]$$

$$=\left[\begin{array}{ccccc}{E}_{11}& \mathbb{O}& \dots & \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \dots & \mathbb{O}& \mathbb{O}\\ \dots & \dots & \dots & \dots & \dots \\ E_{p-1,1}& E_{p-1,2}^{\left(1\right)}& \dots & E_{p-1,p-1}^{(p-2)}& E_{p-1,p}^{(p-2)}\\ E_{p1}& E_{p2}^{\left(1\right)}& \dots & E_{p,p-1}^{(p-2)}& E_{pp}^{(p-2)}\end{array}\right]N_{p,p-1}$$

(59)

and

$$[E_1{E}_2^{\left(1\right)}\dots E_{p-2}^{(p-3)}{E}_{p-1}^{(p-2)}{E}_p^{(p-1)}]$$

$$=[E_1{E}_2^{\left(1\right)}\dots E_{p-2}^{(p-3)}{E}_{p-1}^{(p-2)}{E}_p^{(p-2)}]N_{p,p-1},$$

(60)

where

$$E_{pp}^{(p-1)}=E_{pp}^{(p-2)}-E_{p,p-1}^{(p-2)}{\left(E_{p-1,p-1}^{(p-2)}\right)}^{†}{E}_{p-1,p}^{(p-2)}$$

(61)

and

$$E_p^{(p-1)}=E_p^{(p-2)}-E_{p-1}^{(p-2)}{\left(E_{p-1,p-1}^{(p-2)}\right)}^{†}{E}_{p-1,p}^{(p-2)}.$$

(62)

As a result, on the basis of (59)–(62), equation (10) is reduced to the equation with the the block-lower triangular matrix as follows

$$\widehat{M}\left[\begin{array}{ccccc}{E}_{11}& \mathbb{O}& \dots & \mathbb{O}& \mathbb{O}\\ E_{21}& E_{22}^{\left(1\right)}& \dots & \mathbb{O}& \mathbb{O}\\ \dots & \dots & \dots & \dots & \dots \\ E_{p-1,1}& E_{p-1,2}^{\left(1\right)}& \dots & E_{p-1,p-1}^{(p-2)}& \mathbb{O}\\ E_{p1}& E_{p2}^{\left(1\right)}& \dots & E_{p,p-1}^{(p-2)}& E_{pp}^{(p-1)}\end{array}\right]=\widehat{E},$$

(63)

where

$$\widehat{M}=[M_1{M}_2\dots M_{p-1}{M}_p]$$

and

$$\widehat{E}=[E_1{E}_2^{\left(1\right)}\dots E_{p-1}^{(p-2)}{E}_p^{(p-1)}].$$

4.1.7. Solution of Equation (10)

Using a back substitution, (63) is represented by p equations as follows:

$$\begin{array}{ccc}\hfill M_p{E}_{pp}^{(p-1)}& =& E_p^{(p-1)},\hfill \\ \hfill M_{p-1}{E}_{p-1,p-1}^{(p-2)}+M_p{E}_{p,p-1}^{(p-2)}& =& E_{p-1}^{(p-2)},\hfill \\ \hfill ⋮\\ \hfill M_2{E}_{22}^{\left(1\right)}+\dots +M_{p-1}{E}_{p-1,2}^{\left(1\right)}+M_p{E}_{p2}^{\left(1\right)}& =& E_2^{\left(1\right)},\hfill \\ \hfill M_1{E}_{11}+\dots +M_{p-1}{E}_{p-1,1}+M_p{E}_{p1}& =& E_1.\hfill \end{array}$$

Therefore, the solution of equation (10) is represented by solutions to the problems

$$\underset{M_p}{min}{\parallel E_p^{(p-1)}-M_p{E}_{pp}^{(p-1)}\parallel }^2,$$

(64)

$$\underset{M_{p-1}}{min}{\parallel (E_{p-1}^{(p-2)}-{\tilde{M}}_p{E}_{p,p-1}^{(p-2)})-M_{p-1}{E}_{p-1,p-1}^{(p-2)}\parallel }^2,$$

(65)

$$⋮$$

$$\underset{M_2}{min}{\parallel (E_2^{\left(1\right)}-\sum _{j=3}^p{\tilde{M}}_j{E}_{j2}^{\left(1\right)})-M_2{E}_{22}^{\left(1\right)}\parallel }^2,$$

(66)

$$\underset{M_1}{min}{\parallel (E_1-\sum _{j=2}^p{\tilde{M}}_j{E}_{j1})-M_1{E}_{11}\parallel }^2.$$

(67)

In (65)–(67), ${\tilde{M}}_p,\dots ,{\tilde{M}}_2$ are solutions to the problems in (64)–(66), respectively. On the basis of [26], the minimal Frobenius norm solutions ${\tilde{M}}_p,\dots ,{\tilde{M}}_1$ to (64)–(67) are given by

$${\tilde{M}}_p=E_p^{\left(p-1\right)}{\left(E_{pp}^{\left(p-1\right)}\right)}^{†},$$

(68)

$${\tilde{M}}_{p-1}=(E_{p-1}^{(p-2)}-{\tilde{M}}_p{E}_{p,p-1}^{(p-2)}){\left(E_{p-1,p-1}^{(p-2)}\right)}^{†},$$

(69)

$$⋮$$

$${\tilde{M}}_2=(E_2^{\left(1\right)}-\sum _{j=3}^p{\tilde{M}}_j{E}_{j2}^{\left(1\right)}){\left(E_{22}^{\left(1\right)}\right)}^{†},$$

(70)

$${\tilde{M}}_1=(E_1-\sum _{j=2}^p{\tilde{M}}_j{E}_{j1})E_{11}^{†}.$$

(71)

Thus, the p-th degree filter requires computation of p pseudo-inverse matrices, $E_{11}^{†}\in {\mathbb{R}}^{q_1\times q_1}$, ${\left(E_{22}^{\left(1\right)}\right)}^{†}\in {\mathbb{R}}^{q_2\times q_2}$, $\dots ,$${\left(E_{pp}^{\left(p-1\right)}\right)}^{†}\in {\mathbb{R}}^{q_p\times q_p}$.

The algorithm that follows is based on formulas (51), (52), (55), (56), (61), (62) and (68)–(71).

Algorithm 1: Solution of equation (10)

Input: p, $E_{kj}$ for $k,j=1,\dots ,p$, and $E_k$ for $k=1,\dots ,p$. Output: $M_j$ for $j=1,\dots ,p$.

Note that the algorithm is quite simple. Its output is constructed from the successive computation of matrices $E_{kj}$ and $E_k$.

• In line 2, the algorithm updates matrices $E_{kj}$ and $E_k$ similar to those for the particular cases given by (51), (52), (55), (56), (61), (62).
• In line 5, matrix $M_p$ is calculated as in (68).
• In line 7, matrices $M_{p-1},\dots ,M_1$ are calculated as in (69)–(71).

Importantly, $E_{ss}^{†}$ is the $q_s\times q_s$ matrix with $q_s<q$, therefore, it requires less computation than $q\times q$ pseudo-inverse $E_{vv}^{†}$ in (7).

In the Matlab code below, $m,q,p$ are as above and Evv.mat, Exv.mat are files containing matrices $E_{vv}$ and $E_{xv}$, respectively.

Listing 1. MATLAB Code for Solving Equation (10).

Example 2.

Here, we wish to numerically illustrate the above algorithm andMatlabcode. To this end, in equation (10), we consider matrices $E_{vv}\in {\mathbb{R}}^{q\times q}$ and $E_{xv}\in {\mathbb{R}}^{m\times q}$ whose entries are normally distributed with mean 0 and variance 1. We choose $m=1000$, $q=1000,\dots ,5000$, $p=1,2,4,5$ and $q_j=q/p$, for $j=1,\dots ,p$.

In Table 2 and Figure 1, for $p=1,2,4,5$, diagrams of the time needed by the above algorithm versus different values of dimension q are presented. Recall, the simulations carry out by the Dell Latitude 7400 Business Laptop (manufactured in 2020). The time decreases with the increase in p. In other words, the computational load needed for the numerical realization of the proposed p-th degree filer decreases with the increase in degree p. Note that the proposed p-th degree filter requires computation of p pseudo-inverse matrices and also matrix multiplications determined by the above algorithm. In particular, for $q=5000$ and $p=5,$ the proposed filter is around four times faster than the filter defined by (7).

4.1.8. The Error Associated with the Filter Represented by (4) and (68)–(71)

Let us denote by

$${\epsilon }_p=\underset{M_1,\dots ,M_p}{min}{∥\mathbf{x}-[M_1,\dots ,M_p]\mathbf{v}∥}_{\Omega }^2$$

(72)

the error associated with the filter determined by (4) and (68)–(71). We wish to analyze the error.

Theorem 4.

The error associated with the filter determined by (4) and (68)–(71) is given by

$${\epsilon }_p={\parallel E_{xx}^{1/2}\parallel }^2-\sum _{j=1}^p{∥E_j^{(j-1)}{\left[{\left(E_{jj}^{(j-1)}\right)}^{1/2}\right]}^{†}∥}^2.$$

(73)

Let ${\mathbf{v}}_p$ be a nonzero vector. Then the error decreases if the degree of the proposed filter p increases, i.e.,

$${\epsilon }_p<{\epsilon }_{p-1},\mathrm{for}p=2,3,\dots .$$

(74)

In particular, if ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$ are such that

$$\sum _{j=1}^p{∥E_j^{(j-1)}{\left[{\left(E_{jj}^{(j-1)}\right)}^{1/2}\right]}^{†}∥}^2>{\parallel E_{xy}{E}_{yy}^{1/2}\parallel }^2$$

(75)

then ${\epsilon }_p$ is less than ϵ, the error associated with the filter given by (3),

$${\epsilon }_p<ϵ.$$

(76)

Proof. 

We have

$${\epsilon }_p={\parallel E_{xx}^{1/2}\parallel }^2-\mathrm{tr}\left\{E_{xv}{E}_{vv}^{†}{E}_{vx}\right\},$$

(77)

where similar to the proofs of (2) and (3), $[{\tilde{M}}_1\dots {\tilde{M}}_p]=E_{xv}{E}_{vv}^{†}$. Note that, as before, $E_1=E_1^{\left(0\right)}$ and $E_{11}=E_{11}^{\left(0\right)}$. Therefore, by (68)–(71),

$$\begin{array}{c}\mathrm{tr}\left\{E_{xv}{E}_{vv}^{†}{E}_{vx}\right\}\hfill \\ =\mathrm{tr}\left\{[{\tilde{M}}_1\dots {\tilde{M}}_p]E_{vx}\right\}\hfill \\ =\mathrm{tr}\{{\tilde{M}}_1{E}_1^T+\dots +{\tilde{M}}_p{E}_p^T\}\hfill \\ =\mathrm{tr}\{(E_1{E}_{11}^{†}-[{\tilde{M}}_2{E}_{21}+\dots +{\tilde{M}}_p{E}_{p1}]E_{11}^{†})E_1^T\hfill \\ +{\tilde{M}}_2{E}_2^T+\dots +{\tilde{M}}_p{E}_p^T\}\hfill \\ =\mathrm{tr}\left\{E_1{E}_{11}^{†}{E}_1^T\right\}+\mathrm{tr}\{{\tilde{M}}_2(E_2^T-E_{21}{E}_{11}^{†}{E}_1^T)+\dots \hfill \\ +{\tilde{M}}_p(E_p^T-E_{p1}{E}_{11}^{†}{E}_1^T)\}+\mathrm{tr}\{[E_2^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†}-{\tilde{M}}_3{E}_{32}^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†}-\hfill \\ \dots -{\tilde{M}}_p{E}_{p2}^{\left(1\right)}{\left(E_{22}^{\left(1\right)}\right)}^{†}]{\left(E_2^{\left(1\right)}\right)}^T+\dots +{\tilde{M}}_p(E_p^T-E_{p1}{E}_{11}^{†}{E}_1^T)\}\hfill \\ =\dots \hfill \\ =\parallel E_1{E}_{11}^{1/2}{\parallel }^2+\parallel E_2{E}_{22}^{1/2}{\parallel }^2+\dots +{\parallel E_{p-1}{E}_{p-1,p-1}^{1/2}\parallel }^2\hfill \\ +\mathrm{tr}\left\{E_p^{\left(2\right)}{\left(E_{pp}^{\left(2\right)}\right)}^{†}{\left(E_p^{\left(2\right)}\right)}^T\right\}\hfill \\ =\parallel E_1{E}_{11}^{1/2}{\parallel }^2+\parallel E_2^{\left(1\right)}{\left[{\left(E_{22}^{\left(1\right)}\right)}^{1/2}\right]}^{†}{\parallel }^2+\dots +{\parallel E_p^{\left(2\right)}{\left[{\left(E_{pp}^{\left(2\right)}\right)}^{1/2}\right]}^{†}\parallel }^2.\hfill \end{array}$$

(78)

As a result, (77) and (78) imply (73). In turn, (73) implies (74). Inequalities (23) and (76) follow from (28), (77) and (78). ▪

Remark 4.

Of course, in practice, the condition in (23) can easily be satisfied by a variation of the dimensions of vectors ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$ and/or by their choice. The following (3) illustrates this observation.

Corollary 3.

If ${\mathbf{v}}_1=\mathbf{y}$ and at least one of ${\mathbf{v}}_2,\dots ,{\mathbf{v}}_p$ is not the zero vector then (76) is always true.

Proof. 

The proof follows directly from (28) and (73). ▪

5. Conclusion

In a number of important applications, such as those in biology and medicine, associated random signals are high-dimensional. Therefore, covariance matrices formed from those signals are high-dimensional as well. Large matrices also appear because of the filter structure. As a result, the numerical realization of the filter that processes high-dimensional signals might be very time consuming, and in some cases run out of computer memory. In this paper, we have presented a method of optimal filter determination which targets the case of high-dimensional random signals. In particular, we have developed a procedure that significantly decreases the computational load needed to numerically realize the filter (see (Section 4.1.7) and (2)). The procedure is based on the reduction of the filter equation with large matrices to its equivalent representation by the set of equations with smaller matrices. The associated algorithm has been provided. The algorithm is easy to implement numerically. We have also shown that the filter accuracy is improved with an increase in filter degree, i.e. in the number of filter parameters (see (4)).

The provided results demonstrate that in high-dimensional random signal settings, a number of existing applications may benefit from using the proposed technique. The Matlab code is given in (Section 4.1.7). In [6], the technique provided here is extended to the case of filtering which is optimal with respect to minimization of both matrices $M_1,\dots ,M_p$ and random vectors ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$.