Sequential Inverse Approximation of a Regularized Sample Covariance Matrix
Tomer Lancewicki
stat.CO






Sequential Inverse Approximation of a Regularized Sample Covariance
Matrix
    Tomer Lancewicki

EBay Inc.

625 6th Ave

New York, NY

Email: tlancewicki@ebay.com
    December 30, 2023
=======================================================================================


One of the goals in scaling sequential machine learning methods pertains
to dealing with high-dimensional data spaces. A key related challenge
is that many methods heavily depend on obtaining the inverse covariance
matrix of the data. It is well known that covariance matrix estimation
is problematic when the number of observations is relatively small
compared to the number of variables. A common way to tackle this problem
is through the use of a shrinkage estimator that offers a compromise
between the sample covariance matrix and a well-conditioned matrix,
with the aim of minimizing the mean-squared error. We derived sequential
update rules to approximate the inverse shrinkage estimator of the
covariance matrix. The approach paves the way for improved large-scale
machine learning methods that involve sequential updates. 





§ INTRODUCTION


The covariance matrix of multivariate data is required in many sequential
machine learning and neural-networks (NN) based applications <cit.>,
including speech recognition <cit.>, deep learning architectures
for image processing and computer vision <cit.>,
stochastic fuzzy NN's <cit.>, pricing option contracts in financial
markets <cit.>, adaptive tracking control problems <cit.>,
detection tasks <cit.>, reinforcement learning <cit.>,
and many others.

In settings where data arrives sequentially, the covariance matrix
is required to be updated in an online manner <cit.>.
Techniques such as cross-validation, which attempt to impose regularization,
or model selection are typically not feasible in such settings <cit.>.
Instead, to minimize complexity, it is often assumed that the covariance
matrix is known in advance <cit.> or that it is restricted to
a specific simplified structure, such as a diagonal matrix <cit.>.
Moreover, when the number of observations n is comparable to the
number of variables p the covariance estimation problem becomes
far more challenging. In such scenarios, the sample covariance matrix
is not well-conditioned nor is it necessarily invertible (despite
the fact that those two properties are required for most applications).
When n≤ p, the inversion cannot be computed at all <cit.>.

An extensive body of literature concerning improved estimators in
such situations exists <cit.>. However, in the absence
of a specific knowledge about the structure of the true covariance
matrix, the most successful approach so far has, arguably, been shrinkage
estimation <cit.>. It has been demonstrated in
<cit.> that the largest sample eigenvalues are systematically
biased upward, and the smallest ones downward. This bias is corrected
by pulling down the largest eigenvalues and pushing up the smallest
ones, toward their grand mean.

The optimal solution of the shrinkage estimator is solved analytically,
which is a huge advantage for deep learning architectures, since a
key factor in realizing such architectures is the resource complexity
involved in their training <cit.>. An example
of such deep architecture is the deep spatiotemporal inference network
(DeSTIN) <cit.>. The latter extensively utilizes the quadratic
discriminant analysis (QDA) classifier under the simplified assumption
that the covariance matrices involved in the process are diagonal.
Such assumption is made in order to avoid additional complexity during
the training and inference processes. It is well known that for a
small ratio of training observations n to observation dimensionality
p, the QDA classifier performs poorly, due to highly variable class
conditional sample covariance matrices. In order to improve the classifiers'
performance, regularization is required, with the aim of providing
an appropriate compromise between the bias and variance of the solution.
It have been demonstrated in <cit.>
that the QDA classifier can be improved tremendously using shrinkage
estimators. The sequential approximated inverse of the shrinkage estimator,
derived in this paper, allows us to utilize the shrinkage estimator
in the DeSTIN architecture with relatively negligible additional complexity
to the architecture. In addition, the relatively simple update rules
pave the way to implement the inverse shrinkage estimator on analog
computational circuits, offering the potential for large improvement
in power efficiency <cit.>. 

The rest of this paper is organized as follows: Section 2 presents
the general idea of the shrinkage estimator. In Section 3, we derived
a sequential update for the shrinkage estimator, while in Section
4, the related approximated inverses are derived. In Section 5, we
conduct an experimental study and examine the sequential update rules.

Notations: we denote vectors in lowercase boldface letters and matrices
in uppercase boldface. The transpose operator is denoted as (·)^T.
The trace, the determinant and the Frobenius norm of a matrix are
denoted as Tr(·), |·|
and ‖·‖ _F, respectively. The identity
matrix is denoted as 𝐈, while 𝐞=[1,1,…,1]^T
is a column vector of all ones. For any real matrices 𝐑_1
and 𝐑_2, the inner product is defined as ⟨𝐑_1,𝐑_2⟩ =Tr(𝐑_1^T𝐑_2),
where ⟨𝐑_1,𝐑_1⟩ =‖𝐑_1‖ _F^2
<cit.>. 




§ SHRINKAGE ESTIMATOR FOR COVARIANCE MATRICES


We briefly review a single-target shrinkage estimator by following
<cit.>, which is generally applied to high-dimensional
estimation problems. Let {𝐱_i} _i=1^n
be a sample of independent identically distributed (i.i.d.)
p-dimensional vectors drawn from a density with a mean μ
and covariance matrix Σ. When the number of observations
n is large (i.e., n≫ p), the most common estimator of Σ
is the sample covariance matrix 

    𝐒_n=1/n-1∑_i=1^n(𝐱_i-𝐦_n)(𝐱_i-𝐦_n)^T,

where 𝐦_n is the sample mean, defined as 

    𝐦_n=1/n∑_i=1^n𝐱_i.

Both 𝐒_n and 𝐦_n are unbiased estimators
of Σ and μ, respectively, i.e.,
E{𝐒_n} =Σ and E{𝐦_n} =μ.
The shrinkage estimator Σ̂(λ_n)
is in the form 

    Σ̂(λ_n)=(1-λ_n)𝐒_n+λ_n𝐓_n

where the target 𝐓_n is a restricted estimator of Σ
defined as 

    𝐓_n=Tr(𝐒_n)/p𝐈.

The work in <cit.> proposed to find an estimator Σ̂(λ_n)
which minimizes the mean squared error (MSE) with respect to
λ_n, i.e., 

    λ_On=min_λ_nE{‖Σ̂(λ_n)-Σ‖ _F^2}

and can be given by the distribution-free formula 

    λ_On=E{⟨𝐓_n-𝐒_n,Σ-𝐒_n⟩}/E{‖𝐓_n-𝐒_n‖ _F^2}.

The scalar λ_On is called the oracle shrinkage coefficient,
since its depends on the unknown covariance matrix Σ.
Therefore, λ_On (<ref>) must be estimated.
The latter can be estimated from its sample counterparts as in <cit.>.
We denote this estimator as λ̂_On. 




§ SEQUENTIAL UPDATE OF THE SHRINKAGE ESTIMATOR


We want to know what happens to Σ̂(λ_n)
(<ref>) when we add an observation 𝐱_n+1,
using only the current knowledge of 𝐒_n, 𝐦_n
and n. Setting 𝐝_n+1=𝐱_n+1-𝐦_n
while using <cit.>, we have the following
update rules for 𝐦_n (<ref>) and 𝐒_n
(<ref>) when an observation 𝐱_n+1 is added

    𝐦_n+1=𝐦_n+1/n+1𝐝_n+1


    𝐒_n+1=n-1/n𝐒_n+1/n+1𝐝_n+1𝐝_n+1^T.



Based on 𝐒_n+1 (<ref>), we can write the update
rule for the target 𝐓_n (<ref>) as


    𝐓_n+1=n-1/n𝐓_n+1/(n+1)p‖𝐝_n+1‖ _F^2𝐈

By using 𝐒_n+1 (<ref>) and 𝐓_n+1
(<ref>), the update rule for the shrinkage estimator Σ̂(λ_n)
(<ref>) can be written as 

    Σ̂(λ_n+1)=𝐆_n+𝐅_n

where 𝐆_n and 𝐅_n defined as 

    𝐆_n=n-1/nΣ̂(λ_n)+(1-λ_n+1)1/n+1𝐝_n+1𝐝_n+1^T

and 

    𝐅_n=1/(n+1)pλ_n+1‖𝐝_n+1‖ _F^2𝐈


    +n-1/n(λ_n-λ_n+1)(𝐒_n-𝐓_n),

respectively. Based on the above update rules, we derive the sequential
update rules for the inverse of the shrinkage estimator.




§ SEQUENTIAL UPDATE FOR THE INVERSE OF THE SHRINKAGE ESTIMATOR


In this section, we derived approximated inverses of the shrinkage
estimator which are updated sequentially and do not involve any matrix
inversion. We start, therefore, from the inverse of the sample covariance
matrix 𝐒_n+1 that can be obtained from the current inverse
of 𝐒_n (<ref>) using the Sherman-Morrison-Woodbury
matrix identity <cit.> as 

    𝐒_n+1^-1=n/n-1(𝐒_n^-1-𝐒_n^-1𝐝_n+1𝐝_n+1^T𝐒_n^-1/n^2-1/n+𝐝_n+1^T𝐒_n^-1𝐝_n+1).

The last update rule can be used only if 𝐒_n is invertible.
It will not be invertible for n≤ p. Since the shrinkage estimator
Σ̂(λ_n) (<ref>)
is a regularized version of 𝐒_n (<ref>),
an inverse exists for any n. This inverse of Σ̂(λ_n)
(<ref>) involves two main steps. The first one is to
update the inverse of 𝐆_n (<ref>) from an inverse
of Σ̂(λ_n) (<ref>).
The second is to update the next step inverse of Σ̂(λ_n)
from 𝐅_n (<ref>) and the inverse of 𝐆_n
(<ref>) calculated in the first step. Suppose, for example,
that the exact inverse of Σ̂(λ_n)
(<ref>), denoted as Σ̂^-1(λ_n),
is known. In the same manner as in 𝐒_n+1^-1 (<ref>),
the inverse for 𝐆_n (<ref>) can be calculated
from Σ̂^-1(λ_n)
as 

    [                                                                                     𝐆_n^-1=; n/n-1(Σ̂^-1(λ_n)-Σ̂^-1(λ_n)𝐝_n+1𝐝_n+1^TΣ̂^-1(λ_n)/n^2-1/n(1-λ_n+1)+𝐝_n+1^TΣ̂^-1(λ_n)𝐝_n+1). ]



Using <cit.>, the exact inverse of Σ̂(λ_n+1)
can be calculated from 𝐆_n^-1 (<ref>) and 𝐅_n
(<ref>) with p iterations 

    (𝐆_n^(i+1))^-1=(𝐆_n^(i)+𝐟_i𝐞_i^T)^-1


    =(𝐆_n^(i))^-1-(𝐆_n^(i))^-1𝐟_i𝐞_i^T(𝐆_n^(i))^-1/1+𝐞_i^T(𝐆_n^(i))^-1𝐟_i,  i=1,…,p

where 𝐟_i and 𝐞_i are the i columns
of 𝐅_n (<ref>) and the identity matrix 𝐈,
respectively. The inverse of Σ̂(λ_n+1)
(<ref>) is equal to the output of the last iteration,
i.e., 

    Σ̂^-1(λ_n+1)=(𝐆_n^(p+1))^-1.

In order to avoid the calculation of p iterations, we can use approximations
for Σ̂^-1(λ_n+1)
(<ref>). The inverse approximations of the shrinkage estimator
are discussed in the following section.




 §.§ Inverse Approximations for the Shrinkage Estimator


We consider two approximations for Σ̂^-1(λ_n+1)
(<ref>). The first approximation is defined as 

    Σ̃_1^-1(λ_n+1)=𝐆̃_n^-1-α_n𝐆̃_n^-1𝐅_n𝐆̃_n^-1

where 

    𝐆̃_n^-1=


    n/n-1(Σ̃_1^-1(λ_n)-Σ̃_1^-1(λ_n)𝐝_n+1𝐝_n+1^TΣ̃_1^-1(λ_n)/n^2-1/n(1-λ_n+1)+𝐝_n+1^TΣ̃_1^-1(λ_n)𝐝_n+1).

The matrix 𝐆̃_n^-1 (<ref>) differs
from 𝐆_n^-1 (<ref>) in the fact that it relies
on the approximated inverse Σ̃^-1(λ_n)
(<ref>), instead of the exact inverse Σ̂^-1(λ_n)
(<ref>). A possible motivation to justify the update rule
(<ref>) stems from the mean value theorem as explained
in <cit.>. Another motivation arises from the Neumann
series <cit.> where Σ̂^-1(λ_n+1)
(<ref>) is approximately equal to Σ̃^-1(λ_n+1)
(<ref>) for α=1 and relatively small 𝐅_n.
We define α_n as the value that minimizes the reconstruction
squared error, i.e., 

    α_n=min_α‖(𝐆̃_n^-1-α𝐆̃_n^-1𝐅_n𝐆̃_n^-1)Σ̂(λ_n+1)-𝐈‖ _F^2

and is equal to 

    α_n=Tr(𝐆̃_n^-1𝐅_n𝐆̃_n^-1Σ̂(λ_n+1)(𝐆̃_n^-1Σ̂(λ_n+1)-𝐈))/‖𝐆̃_n^-1𝐅_n𝐆̃_n^-1Σ̂(λ_n+1)‖ _F^2

Additional simplification can be taken by looking at the last term
in 𝐅_n (<ref>). Under the assumption that the
difference λ_n-λ_n+1 is relatively small, we can
write an approximation for 𝐅_n (<ref>) by neglecting
its last term, i.e., 

    𝐅̃_n=1/(n+1)pλ_n+1‖𝐝_n+1‖ _F^2𝐈

This will lead to the second approximation for Σ̂^-1(λ_n+1)
(<ref>), denoted as 

    Σ̃_2^-1(λ_n+1)=𝐆̃'̃_n^-1-α'_n𝐆̃'̃_n^-1𝐅̃_n𝐆̃'̃_n^-1

where 

    [                                                                                         𝐆̃'̃_n^-1=; n/n-1(Σ̃_2^-1(λ_n)-Σ̃_2^-1(λ_n)𝐝_n+1𝐝_n+1^TΣ̃_2^-1(λ_n)/n^2-1/n(1-λ_n+1)+𝐝_n+1^TΣ̃_2^-1(λ_n)𝐝_n+1) ]

and α'_n is calculated by 

    α'_n=(n+1)pTr(𝐆̃'̃_n^-2Σ̂(λ_n+1)(𝐆̃'̃_n^-1Σ̂(λ_n+1)-𝐈))/λ_n+1‖𝐝_n+1‖ _F^2‖𝐆̃'̃_n^-2Σ̂(λ_n+1)‖ _F^2

We examine these two approximations in the following section.




§ EXPERIMENTS


In this section we implement and evaluate the sequential update of
the inverse shrinkage estimator. As in <cit.>, we assume that
the observations are i.i.d Gaussian vectors. In order to study the
estimators performance, an autoregressive covariance matrix Σ
is used. We let Σ be the covariance matrix of
a Gaussian AR(1) process <cit.>, denoted by 

    Σ_AR={σ_ij=r^|i-j|} .

As in <cit.>, we use r=0.5. In all simulations, we
set p=50 and let n range from 1 to 30. Each simulation is repeated
200 times and the average values are plotted as a function of n.
The experimental results are summarized in box plots. On each box,
the central mark is the median, the edges of the box are the 25th
and 75th percentiles, and the whiskers correspond to approximately
+/–2.7σ or 99.3 coverage if the data are normally
distributed. The outliers are plotted individually.

The reconstruction errors of the approximated inverses Σ̃_1^-1(λ_n)
(<ref>) and Σ̃_2^-1(λ_n)
(<ref>) are defined by 

    e_1(n)=1/p‖Σ̃_1^-1(λ_n)Σ̂(λ_n)-𝐈‖ _F^2,

and 

    e_2(n)=1/p‖Σ̃_2^-1(λ_n)Σ̂(λ_n)-𝐈‖ _F^2,

respectively. These reconstruction errors are normalized with p
since it is the squared Frobenius norm of the identity matrix 𝐈.
We examine the approximated inverse Σ̃_1^-1(λ_n)
(<ref>) and Σ̃_2^-1(λ_n)
(<ref>) where λ_n is equal to λ̂_On
<cit.>. The experimental results for the reconstruction errors
e_1(n) (<ref>) and e_2(n)
(<ref>) are summarized in Fig. 1 and Fig. 2, respectively.
The values of e_1(n) (<ref>) converge on average
to zero as the number of observations n increase. In several simulations,
however, the update rule accumulates error and diverges.







The related reconstruction error e_2(n) (<ref>)
is depicted in Fig. 2. The reconstruction error e_2(n)
(<ref>) does not converge to zero due to its relative simplification
involving the use of 𝐅̃_n (<ref>) instead
of 𝐅_n (<ref>). However, the use of 𝐅̃_n
(<ref>) renders Σ̃_2^-1(λ_n)
(<ref>) much more robust to outliers in comparison
to the first estimator Σ̃_1^-1(λ_n)
(<ref>). In that sense, a relatively small and fixed
reconstruction error can be assumed in order to avoid unexpected outliers.




§ CONCLUSIONS


A key challenge in many large-scale sequential machine learning methods
stems from the need to obtain the covariance matrix of the data, which
is unknown in practice and should be estimated. In order to avoid
additional complexity during the modeling process, it is commonly
assumed that the covariance matrix is known in advanced or, alternatively,
that simplified estimators are employed. In Section 3, we derived
a sequential update rule for the shrinkage estimator that offers a
compromise between the sample covariance matrix and a well-conditioned
matrix. The optimal shrinkage coefficient, in the sense of mean-squared
error, is analytically obtained, which is a notable advantage since
a key factor in realizing large-scale architectures is the resource
complexity involved. In Section 4, sequential update rules that approximate
the inverse shrinkage estimator are derived. The experimental results
in Section 5 clearly demonstrates that the reconstruction errors of
the approximated inverses are relatively small. The sequential update
rules that approximate the inverse of the shrinkage estimator provide
a general result that can be utilized in a wide range of sequential
machine learning applications. Therefore, the approach paves the way
for improved large-scale machine learning methods that involve sequential
updates in high-dimensional data spaces. 

IEEEbib


