Recovery of Missing Data in Correlated Smart Grid Datasets

TL;DR

This paper investigates the recovery of missing data in correlated smart grid datasets using matrix completion, demonstrating that Bayesian SVT outperforms standard SVT by effectively exploiting dataset correlations.

Contribution

It introduces a framework for joint recovery of correlated datasets in smart grids and evaluates the performance of BSVT, showing its superiority over SVT in this context.

Findings

BSVT outperforms SVT in data recovery accuracy.

Correlation between datasets enhances recovery performance.

Fundamental limits characterize the potential of joint data recovery.

Abstract

We study the recovery of missing data from multiple smart grid datasets within a matrix completion framework. The datasets contain the electrical magnitudes required for monitoring and control of the electricity distribution system. Each dataset is described by a low rank matrix. Different datasets are correlated as a result of containing measurements of different physical magnitudes generated by the same distribution system. To assess the validity of matrix completion techniques in the recovery of missing data, we characterize the fundamental limits when two correlated datasets are jointly recovered. We then proceed to evaluate the performance of Singular Value Thresholding (SVT) and Bayesian SVT (BSVT) in this setting. We show that BSVT outperforms SVT by simulating the recovery for different correlated datasets. The performance of BSVT displays the tradeoff behaviour described by the fundamental limit, which suggests that BSVT exploits the correlation between the datasets in an efficient manner.

Equations

$${\mathbf{M}}^{(s)}\stackrel{\Delta }{=}[{\mathbf{m}}_1^{(s)},{\mathbf{m}}_2^{(s)},...,{\mathbf{m}}_N^{(s)}],$$

$${\mathbf{m}}_i^{(s)}\stackrel{\Delta }{=}[m_{i,1}^{(s)},m_{i,2}^{(s)},...,m_{i,M}^{(s)}{]}^{\mathsf{T}}\in {\mathbb{R}}^M.$$

$${\mathbf{m}}_i^{(s)}\sim \mathcal{N}(\text{boldmath}{\mu }_s,\text{boldmath}{\Sigma }_s),$$

$${\bf M}\stackrel{{\scriptstyle\Delta}}{{=}}\left[\begin{array}[]{c}{\bf M}_{1}\\ {\bf M}_{2}\end{array}\right].$$

$$\begin{array}{r}\mathbf{M}=[{\mathbf{m}}_1,{\mathbf{m}}_2,...,{\mathbf{m}}_N],\end{array}$$

$$\begin{array}{r}{\mathbf{m}}_i\sim \mathcal{N}(0,\mathbf{\Sigma }).\end{array}$$

$$\hbox{\boldmath$\Sigma$}\stackrel{{\scriptstyle\Delta}}{{=}}\left[\begin{array}[]{c c}\hbox{\boldmath$\Sigma$}_{11}&\psi\hbox{\boldmath$\Sigma$}_{11}\\ \psi\hbox{\boldmath$\Sigma$}_{11}&\hbox{\boldmath$\Sigma$}_{22}\\ \end{array}\right],$$

$$\begin{array}{r}(\text{boldmath}{\Sigma }_{ll}{)}_{i,j}\stackrel{\Delta }{=}{\rho }^{\frac{1}{{\zeta }_{ll}}|i-j|},\end{array}$$

$$\text{boldmath}{\Sigma }_{ll}=\text{Toeplitz}(1,\dots ,{\upsilon }_{ll}),$$

$$\begin{array}{r}{\upsilon }_{ll}={\rho }^{\frac{1}{{\zeta }_{ll}}(M-1)}.\end{array}$$

$${\mathbf{R}}_l={\mathbf{M}}_l+{\mathbf{N}}_l,$$

$$({\mathbf{N}}_l{)}_{i,j}\sim \mathcal{N}(0,{\sigma }_{{\mathbf{N}}_l}^2),$$

$${\Omega }_l\stackrel{\Delta }{=}\{(i,j):({\mathbf{R}}_l{)}_{i,j}\text{ is observed}\}.$$

$$P_{{\Omega }_l}({\mathbf{R}}_l)=\{\begin{array}{ll}({\mathbf{R}}_l{)}_{i,j},& (i,j)\in {\Omega }_l,\\ 0,& \text{otherwise}.\end{array}$$

$$P_{\Omega}({\bf R})=\left[\begin{array}[]{c}P_{\Omega_{1}}({\bf R}_{1})\\ P_{\Omega_{2}}({\bf R}_{2})\end{array}\right]\in\mathbb{R}^{2M\times N},$$

$$\textnormal{NMSE}\left({\bf M};g\right)=\frac{\mathbb{E}\left[\|{\bf M}-g\big{(}P_{\Omega_{1}}({\bf R}_{1}),P_{\Omega_{2}}({\bf R}_{2})\big{)}\|^{2}_{F}\right]}{\|{\bf M}\|_{F}^{2}},$$

$$\underset{\mathbf{X}}{\mathrm{minimize}}$$

$$P_{\Omega }(\mathbf{X})=P_{\Omega }(\mathbf{M}),$$

$$\underset{\mathbf{X}}{\mathrm{minimize}}$$

$$P_{\Omega }(\mathbf{X})=P_{\Omega }(\mathbf{M}),$$

$$\begin{cases}\mathbf{X}^{(k)}=D_{\tau}(\mathbf{Y}^{(k-1)}),\\ \mathbf{Y}^{(k)}=\mathbf{Y}^{(k-1)}+\delta_{s}\big{(}P_{\Omega}(\mathbf{M})-P_{\Omega}(\mathbf{X}^{(k)})\big{)},\\ \end{cases}$$

$$\mathbf{Z}=\mathbf{M}+\mathbf{W},$$

$$(\mathbf{W}{)}_{i,j}\stackrel{iid}{\sim }\mathcal{N}(0,{\sigma }_{\mathbf{Z}}^2),$$

$${\mathbf{Z}}^{(k)}={\mathbf{Y}}^{(k)}+{\mathbf{L}}^{(k)},$$

$$\begin{array}{rl}\text{SURE}(D_{\tau })(\mathbf{Z})=& -2MN{\sigma }_{\mathbf{Z}}^2+\sum _{i=1}^N\text{min}({\tau }^2,{\sigma }_i^2(\mathbf{Z}))\\ & +2{\sigma }_{\mathbf{Z}}^2\text{div}(D_{\tau }(\mathbf{Z})),\end{array}$$

$$\begin{split}\textnormal{div}(D_{\tau}({\bf Z}))=&\sum_{i=1}^{N}\bigg{[}\mathds{1}(\sigma_{i}({\bf Z})>\tau)+(2M-N)\frac{(\sigma_{i}({\bf Z})-\tau)_{+}}{\sigma_{i}({\bf Z})}\bigg{]}\\ &+2\sum_{i\neq j,i,j=1}^{N}\frac{\sigma_{i}({\bf Z})(\sigma_{i}({\bf Z})-\tau)_{+}}{\sigma_{i}^{2}({\bf Z})-\sigma_{j}^{2}({\bf Z})},\end{split}$$

$$\text{max}(r_1,r_2)\le r\le r_1+r_2.$$

$$r=r_1+r_2-d,$$

$$d=\text{dim}[\mathcal{C}({\mathbf{M}}_1^{\mathsf{T}})\cap \mathcal{C}({\mathbf{M}}_2^{\mathsf{T}})],$$

$$0\le d\le \text{min}(r_1,r_2),$$

$$k_1>(M+N-r_1)r_1,$$

Full text

Index Terms— smart grid, matrix completion, missing data recovery, correlated data

1 Introduction

The integration of residential low carbon energy sources such as solar or wind power generates bidirectional power flows that affect the stability of the smart grid [1]. The control strategies need to adapt to the new challenges posed by the additional distributed energy sources. In this context, the monitoring procedures are expected to manage the dynamic and unknown scenarios and to provide timely and accurate data describing the state of the grid. For example, the lack of data quality in power systems contributed towards several large-scale blackouts such as the 2003 U.S.-Canadian blackout [2] and the 2003 Italy blackout [3]. In addition, the integration of the Internet of things into the smart grid will significantly increase the number of datasets [4]. In practical scenarios, state estimation and monitoring systems face challenges like data injection attacks [5], [6], [7], [8], [9] or missing data [10], [11], [12]. Telemetry errors such as sensor failures or communication issues lead to incomplete sets of observations that do not fully describe the state of the grid. Therefore, it is vital to estimate the missing data based on the available observations. For instance, accurate measurements are necessary to implement centralized control schemes for voltage regulation in distribution systems [13].

Matrix completion (MC) is proposed in [14] as technique to recover missing data from partial observations. MC-based recovery exploits the fact that correlated state variable vectors give rise to approximately low rank data matrices. Specifically, in a convex optimization context, a low rank matrix is estimated given that a sufficient fraction of the entries is observed. See for instance [15] and [16]. However, when the number of observations is insufficient, the recovery of the data matrix is not possible. A potential way forward in this case is to attempt a joint recovery of multiple datasets by exploiting the fact that when datasets are correlated the rank of the resulting joint dataset grows in a sub-additive fashion. When the number of observations in one dataset is limited, this approach allows the estimation process to incorporate datasets produced by other sources in the system.

A framework for jointly recovering multiple datasets is provided in [17] where the MC setting is extended to the tensor case. Moreover, the singular value decomposition is extended to the tensor case in [18] which leads to the development of a tensor nuclear norm based algorithm in [19]. Alternatively, a collective MC framework is proposed in [20] to exploit the correlation between matrices with shared structure. However, the common structure constraint does not allow for sufficient generality in the definition of the correlation structure between datasets in a smart grid context.

This paper proposes an estimation setting in which data from multiple datasets is combined into a single data matrix that is recovered using MC-based algorithms. This allows the recovery process to exploit not only correlations within a dataset but also between datasets in the joint estimation paradigm. Specifically, the correlation between datasets is leveraged to facilitate the recovery when the number of observations in one dataset is limited. In addition, the fundamental limit of the joint recovery setting for two correlated datasets is characterized within an MC framework and, based on the geometry dictated by the fundamental limit, the joint recovery performance of two MC-based algorithms is benchmarked for different levels of correlation between the combined datasets. Numerical results show that the recently proposed Bayesian Singular Value Theresholding (BSVT) algorithm [12] is more effective in exploiting the correlation between datasets when compared to the Singular Value Theresholding (SVT) algorithm [21].

2 System Model

Consider an electricity distribution system with $N$ low voltage (LV) feeders. At the head of each feeder, a sensing unit measures various electrical magnitudes, e.g., voltage, intensity, active and reactive power at given time instants. These measures comprise the state variables that the operator uses for control, monitoring, and management purposes. The set of observations available to the operator is incomplete and corrupted by noise. The operator estimates the missing data based on the available observations. In the following, the analysis is carried out for a particular electrical magnitude, i.e., phase voltage.

2.1 Source Model

For a given phase voltage state variable, let $m_{i,j}^{(s)}$ be the corresponding value on phase $s\in\{A,B,C\}$, at feeder $i\in\{1,2,...,N\}$ and time $j\in\{1,2,...,M\}$. The matrix with the measurements for phase $s$, denoted by ${\bf M}^{(s)}\in\mathbb{R}^{M\times N}$, contains the aggregated measurement vectors from all feeders

[TABLE]

where the measurement vectors are given by

[TABLE]

The resulting data matrices ${\bf M}^{(\text{A})}$, ${\bf M}^{(\text{B})}$, ${\bf M}^{(\text{C})}$ contain the voltage measurements on phase A, B and C respectively, at time instants $1,2,...,M$ for all $N$ feeders.

2.2 Real data model

Real data collected as part of the “Low Voltage Network Solutions” project run by Electricity North West Limited (ENWL) [22], is used in the following to model the statistical structure of the random process governing the phase voltage state variables. The dataset contains voltage measurements of phases A, B and C collected from 200 residential secondary substations across North West of England from June 2013 to January 2014. Each substation generates a daily file that contains the voltage measurements on all three phases.

An analysis of the distribution and sample covariance matrix of the phase A voltage measurements in the LV dataset under consideration is presented in [11]. Therein, it is shown that voltage measurements can be modelled as a multivariate Gaussian random process for $i\in\{1,2,...,N\}$. Specifically, we model the voltage measurements as

[TABLE]

and ${{\bf m}_{i}}$ for $i\in\{1,2,...,N\}$, is a sequence of independent and identically distributed random variables. Moreover, it is also shown in [11] that the sample covariance matrix for phase A exhibits a structure that is approximately Toeplitz. In addition, because the voltage data is correlated, the covariance matrix displays a high correlation across feeders and time instants. It is shown in [12] that the singular value decomposition of a $500\times 500$ matrix with phase A voltage measurements has a large condition number [23].

Part of the LV data collected by ENWL is used to construct two complete data matrices ${\bf M}^{(\text{B})}$ and ${\bf M}^{(\text{C})}$ with $M=N=500$ that contain phase B and phase C voltage measurements from the LV grid. The sample covariance matrix of the data matrix ${\bf M}^{(\text{B})}$ is depicted in Fig. 1(a). As expected, and in agreement with the observation in [11], the sample covariance matrix for the phase voltage data exhibits a structure that is approximately Toeplitz. In addition, when the phase B and phase C data matrices are combined into a single data matrix, i.e., ${\bf M}^{(\text{BC})}=[{\bf M}^{(\text{B})},\,{\bf M}^{(\text{C})}]^{\sf T}$, the resulting sample covariance matrix is depicted in Fig. 1(b). Interestingly, the sample covariance for the combined matrix is a block matrix with four elements where each element exhibits a structure that is approximately Toeplitz. Based on this observation, the following section proposes a general model for correlated voltage datasets generated by different phases in smart grid systems.

2.3 Synthetic data model

A mathematical description of the model used to generate two correlated synthetic datasets follows. Let us denote the data matrix for the first dataset by ${\bf M}_{1}\in\mathbb{R}^{M\times N}$ and the data matrix for the second dataset by ${\bf M}_{2}\in\mathbb{R}^{M\times N}$. In this setting, the combined matrix is denoted by ${\bf M}\in\mathbb{R}^{2M\times N}$ given by

[TABLE]

Hence, the combined state variable matrix is defined as

[TABLE]

where each state variable vector ${\bf m}_{i}\in\mathbb{R}^{2M}$ for $i\in\{1,2,...,N\}$ is generated by a multivariate Gaussian process with $\mathbf{0}$ mean and covariance matrix $\Sigma$, i.e.,

[TABLE]

The covariance matrix $\Sigma$ is a block matrix in which block \hbox{\boldmath\Sigma}_{ll} is a Toeplitz matrix describing the covariance matrix of the dataset $l\in\{1,2\}$. The resulting covariance matrix is given by

[TABLE]

where \hbox{\boldmath\Sigma}_{ll}\in\mathbb{R}^{M\times M} and $\psi\in[0,1]$. In this framework, the elements of \hbox{\boldmath\Sigma}_{ll} are defined as

[TABLE]

where (\hbox{\boldmath\Sigma}_{ll})_{i,j} denotes the entry in row $i$ and column $j$ of the matrix \hbox{\boldmath\Sigma}_{ll} with $i\in\{1,2,\ldots,M\}$, $j\in\{1,2,\ldots,M\}$, $\rho\in(0,1)$ and $\zeta_{ll}$ a design parameter. Hence, the matrix \hbox{\boldmath\Sigma}_{ll} is given by

[TABLE]

where $\upsilon_{ll}\in[0,1)$ obeys

[TABLE]

Fig. 2 describes the system model for the joint recovery of two datasets produced in an LV distribution system. In this setting, each phase voltage data matrix fully describes the state of the grid over $M$ time instants and across $N$ feeders. However, in the acquisition process, part of the measurements are lost and the ones that are available are corrupted by noise. The missing data recovery strategy needs to estimate the actual state of the grid for a noisy subset of observations.

2.4 Acquisition

The phase voltage measurements are assumed to be corrupted by additive white Gaussian noise (AWGN) such that for each dataset the resulting observations are given by

[TABLE]

where $l\in\{1,2\}$ denotes the number of datasets and

[TABLE]

where $i\in\{1,2,...,M\}$ and $j\in\{1,2,...,N\}$. Moreover, it is also assumed that only a fraction of the complete set of observations (entries in ${\bf R}_{l}$) are communicated to the operator. Denote by $\Omega_{l}$ the subset of observed entries of the dataset $l$, i.e.,

[TABLE]

Formally, the acquisition process is modelled by the functions $P_{\Omega_{l}}:\mathbb{R}^{M\times N}\to\mathbb{R}^{M\times N}$ with $l\in\{1,2\}$ and

[TABLE]

The observations given by (14) describe all the data from dataset $l$ that is available to the operator for estimation purposes. Therefore, the recovery of the missing data is performed from the observations $P_{\Omega_{l}}({\bf R}_{l})$. As depicted in Fig. 2, the acquisition step is performed independently for each dataset. After the acquisition step, the available observations from each dataset are combined into a single data matrix, i.e.,

[TABLE]

where $\Omega$ denotes the combined set of available observations from the two datasets. The resulting matrix $P_{\Omega}({\bf R})$ is used for estimation purposes in the joint recovery paradigm.

2.5 Estimation

The estimation process for the combined matrix of measurements, based on the available observations from each dataset is modelled by the function $g:\mathbb{R}^{2M\times N}\to\mathbb{R}^{2M\times N}$, where $\Omega$ denotes the combined set of available observations from both datasets. The estimate \widehat{{\bf M}}=g\big{(}P_{\Omega_{1}}({\bf R}_{1}),P_{\Omega_{2}}({\bf R}_{2})\big{)} is obtained by solving an optimization problem based on a given optimality criterion. In the following, the optimality criterion is the normalized mean square error (NMSE) given by

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm.

3 Recovering missing data using matrix completion

Given a matrix ${\bf M}$ of size $2M\times N$, and observations $P_{\Omega}({\bf M})$, the recovery of the missing entries is not feasible in the general case. However, when ${\bf M}$ is low rank or approximately low rank, it is shown in [14] that if the entries on $\Omega$ are sampled uniformly at random, the missing entries are recovered with high probability by solving the following optimization problem:

[TABLE]

where $\lVert\mathbf{X}\rVert_{*}$ denotes the nuclear norm of the matrix ${\bf X}$. To simplify the notation, let us assume that $2M\geq N$. We proceed to present the two MC-based algorithms used to assess the joint recovery performance. Namely, the SVT algorithm proposed in [21] and the BSVT approach presented in [12].

3.1 Singular Value Theresholding

SVT is an MC-based algorithm [21] which produces a sequence of matrices ${\bf X}^{(k)}$ that converges to the unique solution of the following optimization problem:

[TABLE]

Note that when $\tau\to\infty$, the optimization problem in (18) converges to the nuclear norm minimization problem in (17). The iterations of the SVT algorithm are:

[TABLE]

where $\mathbf{Y}^{(0)}=\mathbf{0}$ is used for initialization, $\delta_{s}$ is the step size that obeys $0<\delta_{s}<2$, and the soft-thresholding operator, $D_{\tau}$ that shrinks the singular values of $\mathbf{Y}^{(k-1)}$ towards zero [21].

Interestingly, the choice of $\tau$ is important to guarantee a successful recovery, since large values guarantee a low-rank matrix estimate but for values larger than $\underset{i}{\textnormal{max}}\,(\sigma_{i}({\bf Y}))$ all the singular values vanish. In [21], the proposed threshold is $\tau=5N$. However, simulation results presented in [11] show that $\tau=5N$ gives suboptimal performance when the number of missing entries is large. The main shortcoming of the SVT algorithm is the lack of guidelines for tuning the threshold $\tau$. This problem is addressed in [12] where a new algorithm is proposed to adapt the recovery to the dataset by leveraging knowledge of the second order statistics.

3.2 Bayesian Singular Value Theresholding

BSVT is an MC-based algorithm [12] that is able to optimize the value of $\tau$ at each iteration using additional prior knowledge in the form of second order statistics. The optimization of the soft-theresholding step is performed using Stein's unbiased risk estimate (SURE) [24] for which a closed-form expression is presented in [25]. However, the result therein pertains to input matrices $\mathbf{Z}$ that accept the following model:

[TABLE]

where the entries of ${\bf W}$ are

[TABLE]

where $\sigma_{\mathbf{Z}}^{2}$ is the variance of the $(\mathbf{W})_{i,j}$ entries with $i\in\{1,2,\ldots,2M\}$ and $j\in\{1,2,\ldots,N\}$. Using the prior knowledge in the form of the second order statistics, the BSVT algorithm, computes the matrix ${\bf Z}$ at iteration $k$, i.e., ${\bf Z}^{(k)}$, as

[TABLE]

where ${\bf Y}^{(k)}$ is defined in (19) and ${\bf L}^{(k)}$ is the linear minimum mean square error (LMMSE) estimate. Consequently, the incorporation of the LMMSE step into the structure of the BSVT algorithm facilitates the use of SURE [24] which is given by

[TABLE]

where $\sigma_{i}({\bf Z})$ is the $i$-th singular value of ${\bf Z}$ for $i\in\{1,2,\ldots,N\}$. A closed-form expression for the divergence of this estimator is obtained in [25]. For the case in which ${\bf Z}\in\mathbb{R}^{2M\times N}$, the divergence is given by

[TABLE]

where $\mathds{1}(\cdot)$ denotes the indicator function. In addition, when ${\bf Z}$ has repeated singular values, the divergence in (24) is set to be zero. The proposed algorithm approximates $\sigma_{{\bf Z}}^{2}$ with the weighted sum of the noise in $\Omega$ and in $\Omega^{c}$. A detailed description of the BSVT algorithm is presented in [12] but we reproduce the algorithm below to aid with the presentation. Note that ${D}_{\textnormal{LMMSE}}$ represents the average noise per entry in $\Omega^{c}$. The main advantage of the BSVT algorithm is that the threshold is optimized at each iteration. This is achieved by incorporating the prior knowledge about the matrix in the form of second order statistics via the introduction of the SURE and LMMSE steps. Admittedly, this approach requires additional knowledge that is not necessary when using the SVT algorithm. However, it is shown in [12] that the introduction of the prior knowledge enables a robust recovery of the missing entries.

4 Joint recovery of missing data in two datasets

In this section, an estimation framework for recovering missing data from different datasets is proposed. The estimation framework facilitates considering the correlation between datasets for a wide range of correlation structures. We begin by noting that in the joint recovery case, there are two types of correlation between the entries of the combined matrix. First, the intra-correlation that refers to the correlation between the entries within each dataset. This is the type of correlation that is exploited in the independent recovery scenario, i.e., when the missing entries from each dataset are recovered using only available observations from that dataset. Second, the cross-correlation defined as the correlation between the data points from the two different datasets. In contrast to the independent recovery case, a joint recovery technique needs to account for both types of correlation. By considering the cross-correlation, the recovery process leverages on other types of data in order to recover the datasets with limited available observations.

In an MC setting, the minimum number of observations required depends on the size and the rank of the matrix [14]. The combination of the datasets into a single matrix increases the size of the matrix, and therefore, the fundamental limit for the joint recovery case depends on the tradeoff between the size and the rank of the combined matrix, and the number of observations available for each dataset. Note that the rank of the combined matrix depends on both the intra and the cross-correlation. The following lemma provides lower and upper bounds for the rank of the combined matrix based on the individual rank of the matrices. To that end, let us denote rank of the matrices by $\textnormal{rank}({\bf M}_{1})=r_{1}$, $\textnormal{rank}({\bf M}_{2})=r_{2}$, and $\textnormal{rank}({\bf M})=r$.

Lemma 1**.**

Let ${\bf M}_{1}\in\mathbb{R}^{M\times N}$ and ${\bf M}_{2}\in\mathbb{R}^{M\times N}$. Define the combined matrix {\bf M}=\left[\begin{array}[]{c}{\bf M}_{1}\\ {\bf M}_{2}\end{array}\right]\in\mathbb{R}^{2M\times N}. Then, the following holds:

[TABLE]

Proof.

The rank of the matrix ${\bf M}$ is defined as (3.23(a) in [23])

[TABLE]

where $d$ is the intersection of the row subspaces of the matrices ${\bf M}_{1}$ and ${\bf M}_{2}$ given by

[TABLE]

where ${\cal C}({\bf M})$ denotes the column subspace of the matrix ${\bf M}$. Since the intersection of the row subspaces is bounded by

[TABLE]

the rank of the matrix ${\bf M}$ satisfies (25). ∎

Based on the insight provided by Lemma 1, the intra-correlation determines the rank of the matrices ${\bf M}_{1}$ and ${\bf M}_{2}$ which define the lower and upper bounds on $r$. A smaller value of intra-correlation in one of the datasets results in a larger lower bound for $r$. On the other hand, the cross-correlation governs the value of $r$ within the limits defined by Lemma 1. Indeed, a larger value of cross-correlation results in a value of $r$ that is closer to the lower bound while a smaller value of cross-correlation generates a combined matrix with a rank that is closer to the upper bound. In other words, the intra-correlation defines the limit values of $r$ for which recovery is feasible and the cross-correlation governs the value of $r$ within the limit.

In [26] it is shown that the low rank matrices ${\bf M}_{1}$ and ${\bf M}_{2}$ can be successfully recovered independently when the number of observations for the first matrix, denoted by $k_{1}$, satisfies

[TABLE]

and the number of available observations for the second matrix, denoted by $k_{2}$, obeys

[TABLE]

This result is based on the assumption that for the random matrices ${\bf M}_{1}$ and ${\bf M}_{2}$ there exist the $\sigma$-measures $\mu_{1}$ and $\mu_{2}$, respectively, and that both measures admit a Lebesgue decomposition. For the combined matrix ${\bf M}$, the $\sigma$-measure is obtained as the product of the measures of ${\bf M}_{1}$ and ${\bf M}_{2}$ [27]

[TABLE]

Moreover, since $\mu$ is a $\sigma$-measure it also admits a Lebesgue decomposition [28] and [29]. Hence, the result in [26] applies for the combined matrix ${\bf M}$ without any additional assumptions, i.e.,

[TABLE]

Fig. 3 depicts the inequalities in (29), (30) and (32) that describe the lower bound on the number of observations required to recover the matrices ${\bf M}_{1}$, ${\bf M}_{2}$ and ${\bf M}$, respectively. The bounds divide the $(k_{1},k_{2})$ plane into seven regions that correspond to different recovery scenarios for the independent and joint estimation settings. The seventh region, i.e., ${\cal R}_{7}$, corresponds to the case in which the independent recovery of each dataset is not possible but the joint recovery is feasible. In other words, the existence of region ${\cal R}_{7}$ is equivalent to the case in which it is beneficial to jointly recover the two datasets. This scenario is captured by the following definition.

Definition 1**.**

The joint recovery of two matrices, ${\bf M}_{1}$, ${\bf M}_{2}\in\mathbb{R}^{M\times N}$ of rank $r_{1}$ and $r_{2}$, respectively, is beneficial in the region given by

[TABLE]

Note that when the rank of the combined matrix increases, the line described by $(2M+N-r)r=k_{1}+k_{2}$ is shifted towards larger values which for values of $r$ larger than a given threshold induces an empty region ${\cal R}_{7}$. The value of the threshold is given by the sufficient condition in Theorem 1. It is also worth noting that based on the regions depicted in Fig. 3, the number of cases in which the joint recovery is feasible is larger than the number of cases in which the independent recovery of the datasets is possible. The following theorem provides the necessary and sufficient conditions to guarantee that the joint recovery is beneficial.

Theorem 1**.**

Let ${\bf M}_{1},{\bf M}_{2}\in\mathbb{R}^{M\times N}$, with rank $r_{1}$ and $r_{2}$. Then, the joint recovery of the two matrices requires fewer observations than the independent recovery if

[TABLE]

and the rank of the combined matrix satisfies

[TABLE]

Proof.

The proof for the necessary condition in (34) hinges on the fact that the quadratic inequality on $r$ that describes the cases in which the joint recovery is beneficial, given by

[TABLE]

is satisfied for a value bounded by Lemma 1. This is achieved by showing that the smaller root of the quadratic is contained in the open interval \big{(}\textnormal{max}(r_{1},r_{2}),r_{1}+r_{2}\big{)}.

The sufficient condition in (35) is equivalent to showing that the smaller root of the quadratic in (36) is a strict upper bound for $r$. ∎

Note that the necessary condition from Theorem 1, that is the inequality in (34) depends only on the matrices ${\bf M}_{1}$ and ${\bf M}_{2}$ and does not depend on the combined matrix ${\bf M}$. In contrast, the sufficient condition in (35) provides an upper bound for the rank of the combined matrix such that the total number of observations required for the joint recovery is fewer when compared to the independent recovery case. Theorem 1 provides a necessary and sufficient condition for the joint recovery of two data matrices to be beneficial.

5 Numerical results

This section presents a numerical evaluation of the joint recovery performance for two datasets. The matrices used for the simulations are generated using the model described in Section 2.3 and the size of the matrices ${\bf M}_{1}$ and ${\bf M}_{2}$ is fixed such that $M=50$ and $N=100$, respectively. Hence, the joint matrix ${\bf M}$ is a square matrix of size $100$. The range of rank values selected aims to characterize the joint recovery in two scenarios: when the two combined matrices have the same rank, i.e., $r_{1}=6$ and $r_{2}=6$, and when the ratio between the two rank values is small, i.e., $r_{1}=6$ and $r_{2}=9$. To facilitate this, the synthetic data model presented in Section 2.3 is used to generate correlated data matrices with the rank values of interest.

5.1 Simulation framework

Using the mathematical model defined in Section 2.3, the covariance matrix for the combined matrix ${\bf M}$ is given by (7) where the intra-correlation between the state variables in ${\bf M}_{1}$ is modelled by the $\upsilon_{11}$ parameter in (10), the intra-correlation between the state variables in ${\bf M}_{2}$ is modelled by the $\upsilon_{22}$ parameter in (10) and the cross-correlation between ${\bf M}_{1}$ and ${\bf M}_{2}$ is modelled by $\psi$ in (7). The numerical analysis shows that a larger value of $\upsilon_{ll}$ results in a more correlated matrix ${\bf M}_{l}$ and consequently a smaller value for $r_{l}$, where $l\in\{1,2\}$. Moreover, the cross-correlation between ${\bf M}_{1}$ and ${\bf M}_{2}$ increases with the value of $\psi$ which leads to a smaller value for $r$ within the bounds defined by Lemma 1.

The matrix ${\bf M}$ generated using the model in (7) is not exactly low rank. Instead, it can be well approximated by a low rank matrix. Let us denote by $\widetilde{{\bf M}}(r)$ the low rank approximation of rank $r$ obtained by vanishing the smallest $N-r$ singular values of the matrix ${\bf M}$. In the following, $r$ is defined as the minimum value for which the NMSE between the matrix ${\bf M}$ and the low rank approximation of rank $r$, i.e., $\widetilde{{\bf M}}(r)$, is below $10^{-3}$. Consequently, the model in (7) is used to generate data matrices ${\bf M}$ such that the low rank approximations $\widetilde{{\bf M}}_{1}(r_{1})$, $\widetilde{{\bf M}}_{2}(r_{2})$ and $\widetilde{{\bf M}}(r)$ have the intended ranks. Moreover, the low rank approximation of the combined matrix, i.e., $\widetilde{{\bf M}}(r)$, is used to evaluate the numerical performance of both BSVT and SVT in exploiting the correlation between the two datasets. The recovery performance of both algorithms is averaged over ten realizations of $\Omega$, where the locations of the available entries are sampled uniformly at random in each dataset.

In the following, a numerical analysis for the joint recovery performance of SVT and BSVT is presented for the cases in which the rank values for the combined matrices are: $r_{1}=6$, $r_{2}=6$, $r=9$ and $r_{1}=6$, $r_{2}=9$, $r=10$. The choice of rank for the combined matrix resembles a high cross-correlation case in which the value of $r$ satisfies the condition imposed by Theorem 1. We consider a low noise regime for which SNR=$50$ dB in both datasets to emphasize the impact of the intra and cross-correlation in the recovery process, where the SNR in dataset $l\in\{1,2\}$ is defined as

[TABLE]

where \hbox{\boldmath\Sigma}_{ll} is defined in (9) and $\sigma_{{\bf N}_{l}}^{2}$ is described in (12). A wider range of rank values and noise regimes is presented in [30]. The efficiency in exploiting the cross-correlation between the combined datasets is evaluated by comparing the recovery performance across different sampling regimes in which the number of observations in $\Omega$ is constant but the ratio between the number of available entries in each dataset varies. The cross-correlation is successfully exploited when the recovery error is similar across different sampling regimes.

Fig. 4 depicts the performance of the BSVT algorithm when $r_{1}=6$, $r_{2}=6$, $r=9$ and SNR=$50$ dB. Interestingly, the contour lines for the $10^{-4}$ and $10^{-3}$ recovery error exhibit a similar shape to the line depicted by $(2M+N-r)r=k_{1}+k_{2}$ in Fig. 3. This suggests that the BSVT algorithm successfully exploits the cross-correlation in that region and obtains a similar recovery performance tradeoff when the ratio between $k_{1}$ and $k_{2}$ varies for a fixed value of $k_{1}+k_{2}$. In contrast, the contour lines for $10^{-3}$, $10^{-2}$ and $10^{-1}$ SVT recovery error depicted in Fig. 5 exhibit a similar shape to the region ${\cal R}_{4}$ in Fig. 3 which corresponds to the independent recovery area in which the cross-correlation is not exploited. Based on this observation, it is reasonable to assume that SVT is not effective in exploiting the cross-correlation as the recovery error changes with the ratio between $k_{1}$ and $k_{2}$ for a fixed total number of observations.

Fig. 6 depicts the performance of the BSVT algorithm when $r_{1}=6$, $r_{2}=9$, $r=10$ and SNR=$50$ dB. In line with the case discussed in Fig. 4, the contour lines for the $10^{-4}$ and $10^{-3}$ recovery error exhibit a similar shape to the fundamental limit in Fig. 3. This suggests that the BSVT approach is able to exploit the cross-correlation between the combined datasets in the almost noiseless regime for both rank cases considered. In Fig. 7 the performance of the SVT algorithm is depicted for the case in which $r_{1}=6$, $r_{2}=9$, $r=10$ and SNR=$50$ dB. In this case, the shape of the contour lines for $10^{-3}$ and $10^{-2}$ recovery error is similar to the shape of the region ${\cal R}_{4}$ in Fig. 3 which suggests that the SVT algorithm is not efficient in exploiting cross-correlation. Consequently, the BSVT algorithm is able to exploit the cross-correlation between the combined datasets more effectively when compared to the SVT approach in the almost noiseless regime. The gain in recovery performance is facilitated by the prior knowledge incorporated in the structure of the BSVT algorithm.

6 Conclusion

The fundamental limits for the joint recovery of two datasets have been characterized in terms of the rank of the single and combined data matrices. Theoretical conditions are derived for the case in which the joint recovery of two datasets requires less observations compared to the independent recovery case. Based on the insight provided by the fundamental limit, the number of cases in which the joint recovery is feasible is significantly larger when compared to the independent recovery setting.

A model for correlated datasets is proposed. Numerical results show that the correlation between different types of data is exploited by leveraging the information provided by the dataset with fewer missing entries to enable the recovery of the other dataset. Moreover, in contrast to the SVT algorithm, the performance of the BSVT approach matches the geometry imposed by the fundamental limit which suggests that BSVT is indeed better suited to exploit the correlation between datasets.