PrecoG: an efficient unitary split preconditioner for the transform-domain LMS filter via graph Laplacian regularization

TL;DR

PrecoG introduces a data-driven, graph Laplacian-based transform to enhance the convergence of transform-domain LMS filters, outperforming existing methods by adaptively learning the optimal preconditioning transformation from input data.

Contribution

The paper proposes PrecoG, a novel method that learns a unitary transform via graph Laplacian regularization to improve LMS filter convergence without prior knowledge of input data class.

Findings

PrecoG significantly improves the condition number over state-of-the-art techniques.

PrecoG effectively generalizes as a split preconditioner for linear systems.

The method demonstrates efficacy in Hebb-LMS settings.

Abstract

Transform-domain least mean squares (LMS) adaptive filters encompass the class of algorithms where the input data are subjected to a data-independent unitary transform followed by a power normalization stage as preprocessing steps. Because conventional transformations are not data-dependent, this preconditioning procedure was shown theoretically to improve the convergence of the LMS filter only for certain classes of input data. However, in reality if the class of input data is not known beforehand, it is difficult to decide which transformation to use. Thus, there is a need to devise a learning framework to obtain such a preconditioning transformation using input data prior to applying on the input data. It is hypothesized that the underlying topology of the data affects the selection of the transformation. With the input modeled as a weighted graph that mimics neuronal interactions, PrecoG obtains the desired transform by recursive estimation of the graph Laplacian matrix. Additionally, we show the efficacy of the transform as a generalized split preconditioner on a linear system of equations and in Hebb-LMS settings. In terms of the improvement of the condition number after applying the transformation, PrecoG performs significantly better than the existing state-of-the-art techniques that involve unitary and non-unitary transforms.

Figures

Equations

$$R_N=E[{\mathbf{x}}_k{\mathbf{x}}_k^H]=\left(\begin{array}{cccc}1& \rho & \cdots & {\rho }^{n-1}\\ \rho & 1& \cdots & {\rho }^{n-2}\\ ⋮& ⋮& \ddots & ⋮\\ {\rho }^{n-1}& {\rho }^{n-2}& \cdots & 1\end{array}\right)$$

$$\underset{N⟶\infty }{\lim }\mathrm{rank}(\eta )=0.\text{vspace}\ast -0.2cm$$

$$\underset{N\to \infty }{\lim }\det ({\tilde{S}}_N-\lambda {\mathbb{I}}_N)=0,\text{vspace}\ast -0.2cm$$

$$\underset{N\to \infty }{\lim }\det \left({\Delta }_{{\tilde{R}}_N}^{-1/2}{\tilde{R}}_N{\Delta }_{{\tilde{R}}_N}^{-1/2}-\lambda {\mathbb{I}}_N\right)=0.\text{vspace}\ast -0.2cm$$

$$\underset{N\to \infty }{\lim }\det \left({\tilde{R}}_N-\lambda {\Delta }_{{\tilde{R}}_N}\right)=0,\text{vspace}\ast -0.2cm$$

$$U_N=\underset{U_N\in O(N)}{\mathrm{argmin}}\left[\det \left({\tilde{R}}_N-\lambda {\Delta }_{{\tilde{R}}_N}\right)\right].\text{vspace}\ast -0.2cm$$

$$U_N=\underset{U_N\in O(N)}{\mathrm{argmin}}\Vert {\tilde{R}}_N-\lambda {\Delta }_{{\tilde{R}}_N}{\Vert }_F^2.\text{vspace}\ast -0.2cm$$

$$U_N$$

$$\displaystyle w_{i}^{t+1}=w_{i}^{t}-\mu Tr\Big{(}\big{[}\frac{\partial E(p)}{\partial U_{N}}\big{]}^{T}\frac{\partial U_{N}}{\partial w_{i}}\Big{)};0<\mu<1.$$

$$\displaystyle\frac{\partial u_{k,l}}{\partial w_{i}}~{}=~{}Tr\Big{(}\frac{\partial u_{k,l}}{\partial L}\Theta_{i}\Big{)}~{}=~{}Tr\Big{(}\big{[}\frac{\partial L}{\partial u_{k,l}}\big{]}^{-T}\Theta_{i}\Big{)},$$

$$L=U_N\Gamma U_N^T⟹\frac{\partial L}{\partial u_{kl}}=U_N\Gamma J^{mn}+J^{nm}\Gamma U_N^T.$$

$$E_N(\mathbf{w},{ϵ}_1,{ϵ}_2,\beta )=E(p)+\beta ({\mathbf{w}}^T\mathbf{w}-1).$$

$$\frac{\partial E_N}{\partial \mathbf{w}}=\frac{\partial E(p)}{\partial \mathbf{w}}+2\beta \mathbf{w}$$

$$\displaystyle w_{i}^{t+1}=w_{i}^{t}\big{(}1-2\beta\big{)}-\mu Tr\Big{(}\big{[}\frac{\partial E(p)}{\partial U_{N}}\big{]}^{T}\frac{\partial U_{N}}{\partial w_{i}}\Big{)}.$$

$$\displaystyle\frac{\partial\boldsymbol{u}_{i}}{\partial a_{mn}}=\sum_{j\neq i}\frac{1}{(\lambda_{i}-\lambda_{p})}\big{<}\frac{\partial A}{\partial a_{mn}},\boldsymbol{u}_{p}\big{>}\boldsymbol{u}_{p};~{}\lambda_{i}\neq\lambda_{p}.$$

$$\frac{\partial {\mathbf{u}}_i}{\partial w_j}$$

$$w_i(n+1)$$

$$E(W(n+1))$$

$$\displaystyle E(W(n))\big{(}1-2\mu R_{xx}\big{)}$$

$$v_i(n)$$

$$P_i(n)$$

$$W(n+1)$$

$$e(n)$$

$$y(n)$$

$$M(ϵ)$$

$$\frac{\partial M(ϵ)}{\partial U}$$

$$\frac{\partial E(p)}{\partial U_n}$$

$$\frac{\partial M(ϵ)}{\partial U}$$

$$A(\tau ){\mathbf{u}}_i(\tau )$$

$$\dot{A}(\tau ){\mathbf{u}}_i(\tau )+A(\tau ){\dot{\mathbf{u}}}_i$$

$$\displaystyle\big{<}\dot{A}\boldsymbol{u}_{i},\boldsymbol{u}_{i}\big{>}+\big{<}A\dot{\boldsymbol{u}}_{i},\boldsymbol{u}_{i}\big{>}$$

Full text

I Introduction

In 1960, Bernard Widrow and Ted Hoff [1] proposed a class of least mean squares (LMS) algorithms to recursively compute the coefficients of an N-tap finite impulse response (FIR) filter that minimizes the output error signal. This computation is achieved by a stochastic gradient descent approach where the filter coefficients are evaluated as a function of the current error at the output. The LMS algorithm and its variants were subsequently used in myriads of applications, including echo cancellation [2, 3], inverse modeling [4], system identification, signal filtering [5, 6] and several others.

Two of the major issues with this approach are the convergence speed and stability. The filter coefficients (or weights) converge in mean while showing small fluctuation in magnitude around the optimal value. The convergence speed depends on the condition number of the autocorrelation matrix of the input, where a condition number close to unity connotes a fast and stable convergence. Later, adaptive algorithms, such as LSL (least square lattice) and GAL (gradient adaptive lattice) [7, 8] filters were designed to achieve faster convergence, immunity to poor condition number of input autocorrelation matrix, and better finite precision implementation compared to the LMS filter. However, these stochastic gradient filters may sometimes produce significant numerical errors, and the convergence is poor compared to recursive least squares (RLS) filters [9]. Due to the fact that the nature of the autocorrelation matrix is data-dependent, improving its condition number by using a transformation is a way to circumvent the convergence issue in the case of real time data.

In order to obtain well-conditioned autocorrelation matrix of any real world input data, we transform the input a priori, which is popularly known as transform-domain LMS (TDLMS) (See Appendix for the LMS and TDLMS algorithms). The discrete Fourier transform (DFT), discrete cosine transform (DCT) [10, 11, 12, 13]

and others act as suitable off-the-shelf transformations of the input data for such problems. The aforementioned step is immediately followed by a power normalization stage [14, 15] and then used as input to the LMS filter. As a geometrical interpretation, the unitary transformation rotates the mean square error (MSE) hyperellipsoid without changing its shape on the axes of LMS filter weights [15]. The rotation tries to align the axes of the hyperellipsoid to the axes of weights. The power normalization is crucial in enhancing the speed of convergence of LMS filter. The normalization forces the hyperellipsoid to cross all the axes at equal distance from the center of the hyperellipsoid. For a perfect alignment after the transformation, the normalization step turns the MSE hyperellipsoid into a hypersphere [15]. It is important to note that these off-the-shelf transformations are still not data-dependent. They work relatively well for certain classes of data having autocorrelation matrices with special structure (eg. Toeplitz).

The TDLMS filter is flexible as it does not attempt to change the working principles and the architecture of LMS filter. Therefore, the transform-domain module can precede other algorithms, such as RLS, GAL and LSL. Notice that the conventional unitary transformations are independent of the underlying data, hence not optimal in regularizing condition numbers of the autocorrelation matrices of arbitrary real-time datasets. As an example, the DCT has been shown to be near-optimal for Toeplitz matrices. However, the DCT loses its near-optimality in conditioning sparse linear systems.

From a different perspective, the transformation of a matrix, such as autocorrelation matrix to improve the condition number is regarded as a subproblem of preconditioning of matrices [16, 17, 18, 19]. Jacobi [20, 21], Gauss-Seidel [22], approximate inverse [23, 24], incomplete LU factorization [25] preconditioners are examples of such data-dependent transformations that utilize a decomposition of the input matrix consisting of coefficients of linear equations. These algorithms are well-suited for solving linear system of equations. In short, the transformation that improves the performance of the LMS filter can also be applied to solving linear systems of equations.

However, not all the strategies for solving linear systems of equations using matrix-preconditioner are suitable for TDLMS. With careful attention, it can be seen that there is a difference between TDLMS and linear systems in terms of the usage of a preconditioner, as explained below. The preconditioning action is implicit in TDLMS filters (split preconditioner). The autocorrelation matrix is not explicitly used in the LMS architecture. Instead, it is the input data or the transformed input data (in case, we transform the data) that are channeled through the LMS lattice and the update of weights is based on error-correcting learning. While expressing the mean square error at the output as a function of filter weights, it appears that the transformation matrix at the input is attempting to condition the input autocorrelation matrix and the convergence of the LMS algorithm depends on the input autocorrelation matrix. On the other hand, in case of solving linear system of equations, $(Ax=b)$ type, the use of a preconditioner is explicit, which is $M^{-1}Ax=M^{-1}b$, with M as a preconditioner matrix, such as the Gauss-Seidel type.

Let $A$ be a matrix to be preconditioned by another matrix $\zeta$. $\zeta$ is said to be a left, right, and split preconditioner if $\zeta^{-1}A$, $A\zeta^{-1}$, and $U_{1}^{-1}AU_{2}^{-T}$ with $\zeta=U_{1}U_{2}^{T}$ respectively provide improved condition numbers compared to that of A. Gauss-Seidel, incomplete LU, and approximate inverse are examples of a left preconditioner. The transformations in TDLMS algorithm are the unitary split preconditioner type. By unitary, we have \big{(}U_{1}U_{2}^{T}=I\big{)} and $U_{1}=U_{2}=U$. Therefore, it is the unitary split preconditioner matrix that can be applied to TDLMS as well as to solving linear systems of equations. In PrecoG, we aim to learn such unitary split preconditioner from input data. In addition, as the transformation is unitary, the energy of input data remains unchanged after transformation.

Unlike the off-the-shelf, data-independent transformations in TDLMS, the derivation of our unitary split conditioner is motivated by the topology of the structured input data. The topology determines neighborhood relationship between data points, which can be represented using graph-theoretic tools [26, 27, 28]. In recent years, manifold processing and regularization have shown promise in different areas of research [29, 30]. Based on such evidence, we hypothesize that the intrinsic topology of the input data affects the construction of a suitable preconditioner matrix. We estimate a data manifold that provides an alternate set of basis acting as a split preconditioning matrix. The data when projected onto the basis are expected to be decorrelated.

The main contributions of this work are as follows.

• •

PrecoG provides an $optimization$ framework that finds the desired unitary transformation for the preconditioning matrix. We iteratively estimate the underlying topology leveraging graph theory, followed by the computation of desired unitary transformation by using the graph Laplacian and first order perturbation theory. PrecoG significantly outperforms all conventional off-the-shelf transformations.

• •

We show that our approach is equally applicable in preconditioning arbitrary linear systems apart from ameliorating the convergence of LMS filters.

• •

Another advantage of PrecoG is that it can be applied without having a prior knowledge about the process that generates the input data. Estimated autocorrelation matrix serves our purpose.

• •

PrecoG shows its efficacy in supervised and unsupervised learning tasks. Using analogy from Neuroscience, PrecoG is shown to work in the Hebb-LMS paradigm.

To our knowledge, no mathematical framework has been presented so far that could potentially generate such matrices from the data.

I-A Why a graph theoretic approach?

The relevance of using graphs is rooted in neuroscience. The LMS algorithm has been pivotal in the overall error-correcting learning paradigm, mostly because of its mathematical elegance and simple interpretation. Its inventors, Widrow and Hoff, named this adaptive filter as Adaline in 1957. Frank Rosenblatt [31] successfully devised a prototypical version of a supervised binary classifier, called the Perceptron, that superficially mimics the action-potential based firing of a single neuron, where the weight update step utilizes the LMS algorithm. Here, the weight is synonymous with the synaptic strength in neuroscience. That analogy is backed up by the facty that the weights are continuously updated real numbers and the final weight values are ‘endured’ (converged), they have a weak notion of plasticity at the synapse level which is observed when the network is sufficiently trained. In reality, the plasticity in neurons in our brain is a consequence of immensely intricate molecular process [32, 33], and the plasticity of the adaptive filter weights are extremely simple parallel. The Perceptron was extended to multiple layers to effectively tackle problems having nonlinear decision boundaries (such as XOR function). Subsequently, numerous (if not the majority of) neural network models including convolutional networks have used the LMS algorithm to update the network weights.

Tuning synaptic strengths (filter weights) is useful but of limited scope, as it delivers an incomplete portrait of the actual task - adapting the filter weights with the datastream in this work. Researchers have gathered evidence supporting the conclusion that neuronal level dynamics orchestrate complex behavioral functions [34, 35, 36, 37]. In general, neurons inhibit or excite other neurons through specific neuromodulators spewed at the synaptic clefts. In an LMS filter, the filter taps can be conceived as neurons. For example, in our LMS filter setting, we have $N$ neurons that conduct signals, which are weighted (tap weights), summed and passed to another neuron. However, contrary to the neurons in our brain, these filter taps are independent of each other in the traditional LMS implementation. PrecoG generalizes the filter by incorporating interactions among the taps. This interaction is a function of the input data.

The question now becomes: how do these filter taps (equivalently neurons) interact (excite/inhibit) with each other over time in order to expedite the convergence of tap weights (synaptic strengths)? The graph is an efficient tool to encode such interaction [38]. It is worth mentioning a couple of points in this aspect - 1) the dynamics of interaction and the synaptic strength influence each other every time the graph parameters are estimated; 2) we are simply looking for the excitatory/inhibitory drive of each filter tap onto other taps. The resultant graph is simple and undirected. However, a neuron can impart an excitatory or inhibitory drive onto itself (Figure). In biological neurons, these synapses of self-stimulation are called autapses [39]. These are abundant among the fast-spiking inhibitory neurons in the neocortex.

II Graph Theory (in brief)

A graph can be compactly represented by a triplet \big{(}\mathcal{V},\mathcal{E},\boldsymbol{w}\big{)}, where $\mathcal{V}$ is the set of vertices, $\mathcal{E}$ the set of edges, and $\boldsymbol{w}$ the weights of the edges. For a finite graph, $|\mathcal{V}|~{}=~{}N$ which is a finite positive integer, and $|\bullet|$ is the cardinality of a set. By denoting $w_{ij}\in\boldsymbol{w}$ as a real positive weight between two vertices $i$ and $j$ with $i,j\in\{1,2,...,N\}$, the adjacency matrix, $A$ of $\mathcal{G}$ can be given by $a_{ij}=w_{ij}$ with $a_{ii}=0$ for a graph with no self-loop. $A\in\mathcal{R}^{N\times N}$ is symmetric for an undirected graph and can be sparse based on the number of edges. The incidence matrix, $B\in\mathcal{R}^{N\times|\mathcal{E}|}$ of $\mathcal{G}$ is defined as $b_{ij}=1~{}or~{}-1$ where the edge $j$ is incident to or emergent from the vertex $i$. Otherwise $b_{ij}=0$. The graph Laplacian $L\in\mathcal{R}^{N\times N}$, which is a symmetric positive-semidefinite matrix, can be given by $L=BWB^{T}$, where $W$ is a diagonal matrix containing $\boldsymbol{w}$. Determining the topology of input data refers to the estimation of $A$ or $L$ depending on the formulation of the problem at hand.

III Problem statement

Let $x_{k}$ = $[x(k)~{}x(k-1)~{}\dots~{}x(k-N+1)]$ be a $N$-length real valued tap-delayed input signal vector at $k^{th}$ instant. The vector representation is convenient for estimating the input autocorrelation as an ergodic process. Let the autocorrelation matrix, denoted by $R_{N}$ be defined as R_{N}=E\big{(}x_{N}x_{N}^{T}\big{)}. We assume that $\Delta_{Y}$ is the main diagonal of a square matrix $Y$ $(\Delta_{Y}=diag(Y))$. Following this notation, after power normalization the autocorrelation matrix becomes $S_{N}=\Delta_{R_{N}}^{-\frac{1}{2}}R_{N}\Delta_{R_{N}}^{-\frac{1}{2}}$. In general, the condition number of $S_{N}$, $\chi_{S_{N}}$, happens to be significantly large in practical datasets. For example, the condition number of the autocorrelation matrix of a Markov process with signal correlation factor as 0.95 has a $\chi$ of $\mathcal{O}(10^{3})$. Notice that we seek $U_{N}$ to minimize the condition number of $S_{N}$. Let a unitary transformation be $U_{N}$ $(U_{N}U_{N}^{T}=I)$ such that the transformed autocorrelation matrix becomes \tilde{R}_{N}=E\big{[}U_{N}^{T}x_{k}x_{k}^{T}U_{N}\big{]}. Next, $\tilde{R}_{N}$ is subjected to a power normalization stage that produces $\tilde{S}_{N}=\Delta_{\tilde{R}_{N}}^{-\frac{1}{2}}\tilde{R}_{N}\Delta_{\tilde{R}_{N}}^{-\frac{1}{2}}$. Precisely, we want the eigenvalues of \lim_{N\to\infty}\tilde{S}_{N}\in\big{[}1-\epsilon_{2},~{}1+\epsilon_{1}\big{]}, where $\epsilon_{1}$ and $\epsilon_{2}$ are arbitrary constants such that $\chi_{max}\simeq\frac{1+\epsilon_{1}}{1-\epsilon_{2}}$. A schematic of our algorithm is given in Fig. 1.

Let us take an example of a $1^{st}$ order Markov input with the signal correlation factor $\rho$ and autocorrelation matrix, $R_{N}$ as

[TABLE]

It is shown in [15] that $\chi_{S_{N}}\simeq(\frac{1+\rho}{1-\rho})^{2}$, which suggests that $\epsilon_{1}=\rho^{2}+2\rho$, and $\epsilon_{2}=2\rho-\rho^{2}$. After applying the DFT, the condition number becomes $\lim_{N\to\infty}\chi_{\tilde{S}_{N}}=(\frac{1+\rho}{1-\rho})$, which indicates that $\epsilon_{1}=\epsilon_{2}=\rho$. On applying the DCT, $\lim_{N\to\infty}\chi_{\tilde{S}_{N}}=1+\rho$ with $\epsilon_{1}=\rho$ and $\epsilon_{2}=0$.

IV Methodology

The search for $U_{N}$ is carried out through an iterative optimization of an associated cost function. It can be argued from section III that the optimal convergence properties are obtained when $\tilde{S}_{N}$ converges to the identity matrix in the rank zero perturbation sense [15]: $A$ and $B$ with $\eta=A-B$ have the same asymptotic eigenvalue distribution if

[TABLE]

In our case, with $\lambda$ as an eigenvalue, this translates to

[TABLE]

which can be expanded as,

[TABLE]

Eq. (4) can be rearranged as,

[TABLE]

which is a quadratic polynomial of $U_{N}$ as $\tilde{R}_{N}=U_{N}^{T}R_{N}U_{N}$. Given the orthonormality of the eigenvectors $U_{n}$, we can rewrite (5) as

[TABLE]

Here, $O(N)$ is the set of unitary matrices. However, in presence of the determinant in (6), obtaining a closed-form expression for $U_{N}$ is difficult to obtain. To overcome this obstacle, we apply the Frobenius norm in (7).

[TABLE]

The Frobenius norm imposes a stronger constraint compared to (3). In fact, while (5) can be solved if at least one column of $\tilde{R}_{N}-\lambda\Delta_{\tilde{R}_{N}}$ can be expressed as a linear combination of rest of the columns, (7) becomes zero only when $\tilde{R}_{N}-\lambda\Delta_{\tilde{R}_{N}}$ is a zero matrix. In effect, it reduces the search space of $U_{N}$. It is due to the fact that the set of $U_{N}$ that solves eq. (7) is a subset of the $U_{N}$ that are also the solutions of eq. (6).

Here, we address two aspects of the problem. First, (7) attempts to minimize the difference between $\tilde{R}_{N}$, which is $U_{N}^{T}R_{N}U_{N}$, and the scaled diagonal matrix of $\tilde{R}_{N}$. This is necessary because it accounts for the spectral leakage [15] as mentioned later in this section. The second aspect is that (7) seems to diagonalize $\tilde{R}_{N}$ apart from attempting to make the eigenvalues unity only. Here, a solution may be hard to obtain in practice. So, we relax the unity constraint by forcing the eigenvalues to lie within a range $[1-\epsilon_{2},~{}1+\epsilon_{1}]$. By enforcing the constraint, (7) with parameters $p~{}=~{}(\boldsymbol{w},\epsilon_{1},\epsilon_{2})$ becomes

[TABLE]

where $s_{+}=1+\epsilon_{1}$ and $s_{-}=1-\epsilon_{2}$ are the upper and lower bounds respectively for the eigenvalues of $S_{N}$. Using the Hadamard product notation, we can express $\Delta_{\tilde{R}_{N}}$ as $U_{N}^{T}R_{N}U_{N}\circ\mathbb{I}_{N}$. The first two constraints in (8) provide a valley in the space spanned by the eigenvectors in $U_{N}$ if $R_{N}$ is positive definite. The valley exists between two surfaces $(1-\epsilon_{2})U_{N}^{T}R_{N}U_{N}\circ I$ and $(1+\epsilon_{1})U_{N}^{T}R_{N}U_{N}\circ I$. At this point, there might be infinitely-many possible solutions. We add a regularizer on $\boldsymbol{w}$ to obtain an acceptable set of solutions.

The undesired result of the imposed restriction given above is that the convergence time may be significant, and the set of solutions of (8) in terms of $U_{N}$ is significantly smaller than that of (7). Upon approaching a minimum of (7), the speed of gradient descent algorithm drops significantly. Although it is theoretically expected that $\lim_{N\to\infty}\tilde{S}_{N}\in[1-\epsilon_{2},~{}1+\epsilon_{1}]$, in practice, it is difficult to guarantee after a prescribed number of iterative steps.

V Laplacian parametrization

The problem of finding a sub-optimal transform by optimizing eq. (8) is solved by leveraging the graph framework. In this framework, the input data is modeled with a finite, single-connected, simple, and undirected graph endowed with a set of vertices, edges and edge weights. For example, for an LMS filter with $N$ taps, an input signal vector $x_{k}$ has length $N$, which can be represented with $N$ vertices. Basically, each vertex corresponds to one tap of the LMS filter. Using the graph, the unknowns of the optimization in (8) are the number of edges and the associated edge weights. A fully-connected graph with N vertices contains $\frac{N(N-1)}{2}$ edges. A deleted edge can be represented with zero edge weight. We denote the the set of unknown parameters as $\boldsymbol{w}$, which is the set of nonzero weights of the graph.

To find the desired transformation $U_{N}$, the algorithm initializes the weights $\boldsymbol{w}$ with random numbers sampled from a Gaussian distribution with zero mean and unit variance. Let $W$ is the diagonal matrix containing $\boldsymbol{w}$. Then by definition, the graph Laplacian, which is symmetric and positive semidefinite by construction, is given by $L=BWB^{T}$. $B$ is the incidence matrix as mentioned in section II. The spectral decomposition of $L$ provides the matrix of eigenvectors $U$. Finally, $U^{T}x_{k}$ is the transformation that is expected to decorrelate the dataset, which may not be possible due to random initialization. Then the cost function (8) helps update the weights and the search for the desired transformation continues in an iterative fashion until the objective conditions are met.

The cost function in eq. (8) is nonconvex. Therefore, the solution is not guaranteed to be a global optimum. In our work, the required solution is obtained through gradient descent with $\mu$ as the step size parameter. From section II, we obtain that $L=BWB^{T}$. Let $\Theta_{i}=\frac{\partial L}{\partial w_{i}}$, which can be evaluated as $\frac{\partial L}{\partial w_{i}}=B\frac{\partial W}{\partial w_{i}}B^{T}$. Using $\mu$ and $\Theta_{i}$, the update equation is given by

[TABLE]

Here, the computation of $\frac{\partial U_{N}}{\partial w_{i}}$ is perfromed by,

[TABLE]

where by using $J^{mn}=\delta_{mk}\delta_{nl}$, $\frac{\partial L}{\partial u_{kl}}$ can be given by,

[TABLE]

In eq. (9), $t$ is the iteration index, and the computation of $\frac{\partial E(p)}{\partial U_{N}}$ is given in the Appendix. To prevent each $w_{i}$ from erratic values during iteration, we impose $2-$norm on the weight vector $\boldsymbol{w}$. On adding the regularization term to eq. (8), the new cost function becomes

[TABLE]

On differentiating $E_{N}$ with respect to $\boldsymbol{w}$, we get

[TABLE]

Following eq. (13), the iterative update of each weight can be given by,

[TABLE]

However, it is difficult to compute the eq. 9 because $\frac{\partial L}{\partial u_{kl}}$ is non-invertible. The proof is given in Appendix. One can use diagonal compensation by adding $\alpha I(\alpha<0.001)$ to $\frac{\partial L}{\partial u_{kl}}$. However, it produces erroneous solution.

We approach the problem of computing $\frac{\partial U_{N}}{\partial w_{i}}$ using first order eigenvector perturbation. We use the result (see Appendix) that if A is a positive-(semi)definite, symmetric matrix and analytic with respect to its entries, then

[TABLE]

Here, $\boldsymbol{u}_{i}$ and $\lambda_{i}$ are the $i^{th}$ eigenvector and eigenvalue of $A$ respectively. By definition, the graph Laplacian is symmetric and positive semidefinite. Using eq. 15, it can be derived that

[TABLE]

The above formulation can be plugged in \frac{\partial U_{N}}{\partial w_{i}}=\big{[}\frac{\partial\boldsymbol{u}_{1}}{\partial w_{j}}~{}\frac{\partial\boldsymbol{u}_{2}}{\partial w_{j}}~{}...~{}\frac{\partial\boldsymbol{u}_{N}}{\partial w_{j}}\big{]} to be used in eq. 14. In eq. 16, it is evident that the space of incremental changes in an eigenvector with respect to edge weight $w_{i}$ is a spanned by the eigenvectors $\boldsymbol{u}_{i}$.

The data-dependent transformation matrix, $U$ utilizes the autocorrelation matrix of the input data ($\tilde{R}_{N}$ inside $E(p)$). A pertinent question is: how do we get the true autocorrelation matrix for streaming real time data? In response to that question, we want to remind the reader that with the data being channeled to the filter weights, we simultaneously estimate the input autocorrelation matrix using a data window. This estimated matrix is used in deriving the data-dependent transformation, which is applied on the subsequent data. We show in experiments that this procedure significantly improves convergence.

VI Results

We split the result section in two parts. The first part considers systems where the matrix that is to be conditioned is known a priori. The second part presents the performance of PrecoG on simulated on-line data, where autocorrelation matrices are estimated using a window.

Part-I (Linear systems)

We show the effectiveness of our approach in preconditioning different matrices against the preconditioners - DCT, DFT, Jacobi (tridiagonal matrix type), GS ($Gauss-Seidel$) and incomplete LU factorization. To represent the strength of an individual algorithm, we incorporate the condition number of each unconditional matrix with the aforementioned methods. In order to scale the condition numbers obtained from several methods with respect to ours, we define a metric, condition ratio = $\frac{\text{condition number obtained from a method}}{\text{condition number obtained by PrecoG}}$. In some of the results, we compute $log_{10}(\text{condition ratio})$ to mitigate the enormous variance present in the condition ratio scores.

First, we apply our method to precondition a Hilbert matrix [40] which is severely ill-conditioned. Hilbert matrix, $H$ is defined as $H(i,j)=\frac{1}{i+j-1}$. In the experiment, we add a regularizer using \big{(}\alpha\mathcal{I}\big{)} with $0<\alpha\leq 1$ as the regularization coefficient. The condition ratios of the existing algorithms including PrecoG on preconditioning the Hilbert matrices, which are regularized by changing the $\alpha$, are shown in Fig. 2(a). Notice that the $X$-axis is given in $-log_{10}$ scale. Therefore, smaller values at $X$ coordinate indicates higher regularization of the Hilbert matrix. On decreasing the value of $\alpha$, the Hilbert matrix becomes severely ill-conditioned, and the performance of the competitive algorithms except $Gauss-Seidel$ exhibit inconsistent behavior. The DCT performs better near $\alpha=1~{}(-log10(\alpha)=0)$ because of the diagonally dominant nature of the matrix. PrecoG outperformed all the comparative methods.

We also evaluate our algorithm on five different random positive definite matrices with the values taken from a zero-mean and unit-variance Gaussian process. We regularize the matrices to ensure positive-definiteness. It is evident from Fig. 2(b), the condition ratios obtained by applying PrecoG outperformed the DCT, Gauss-Seidel, DFT, and Jacobi transformations.

The condition ratios (in $log10$ scale) with respect to PrecoG on sparse systems of equations are shown in Fig. 2(c). For PrecoG, $log10$ of the condition ratio is zero. So, it is not shown in the plot. The three sparse matrices are random by construction with sparsity \big{(}\frac{\text{number of nonzero elements}}{\text{total number of elements}}\big{)} levels as $[\frac{1}{4},\frac{2}{7},\frac{1}{5}]$ respectively. PrecoG significantly outperformed the conventional transformations.

Part-II (Simulated process)

PrecoG is evaluated on two simulated datasets from two different processes - $1^{st}$ order Markov process and $2^{nd}$ order autoregressive process. The LMS algorithm is supervised on these datasets in a sense that the desired output is known beforehand for each dataset. The data is channeled to the filter in on-line mode, and PrecoG is blinded against the customization of the data. The autocorrelation matrices of both processes are Toeplitz. So far, the DCT is known as the near-optimal preconditioner for $1^{st}$ and $2^{nd}$ order Markov processes. We have also considered the performance of PrecoG on Hebbian-LMS learning [41]. This is an unsupervised form of learning, which is claimed to emulate neuronal learning paradigm.

VI-A $1^{st}$* order Markov or AR(1) process*

We simulate a first order autoregressive process with a set of signal correlation factors,$\rho$. As mentioned in the third claim of our main contributions, PrecoG is shown to performed significantly well without the prior knowledge of asymptotic autocorrelation matrix of the input data.

For each $\rho$, we have simulated $2000$ samples of 1D data and convolved the data with a filter defined by coefficients $h=[1~{}-0.8~{}6~{}3]$ to generate the desired output. The data power is taken as unity. The output is added with a white Gaussian noise of signal power $0.01$. The goal is to converge the filter weight to the impulse function of the convolution filter. Conventional transformations, such as the DCT and the DST are applied at the beginning. However, due to data-dependent nature of PrecoG, we consider a part of initial data to estimate the time averaged autocorrelation function. The final unitary preconditioner of PrecoG depends on the length of the initial data. The step length of the LMS algorithm is set as $0.002$. For a fixed length of initial data, we have created a search space of L2 coefficients and PrecoG learning rate, to find out the transformation.

It is evident from Fig 3(a)that the number of taps in the LMS filter is $4$. We have investigated the relationship between these $4$ taps (neurons) by varying the initial data window. Here, $R-est~{}25$ means that the estimated autocorrelation is averaged over the first $25$ data and the ‘colored’ matrix is the transformation. The color describes the strength of excitatory (red) and inhibitory (blue) drives. The relationship (matrix) is asymmetric. As an example to interpret this asymmetricity, it can be seen that for $R-est~{}25$, neuron 1 inhibits neuron 2, but neuron 2 applies excitatory postsynaptic potential to neuron 1. The diagonal entries of each transformation are non-empty, indicating the presence of autapses. Fig. 3(b) exhibits the transformation of the DCT, the DST and PrecoG, where the PrecoG transformation is obtained from a different instance of AR(1) process with same $h$.

A comparison between the DCT and PreoG in terms of the convergence of filter weights is given in Fig. 4. The DCT is applied to the input data prior to resuming the LMS filter operation. PrecoG functions in two stages. At first, PrecoG computes the transformation using an initial data window (in this example, window length is $200$), during which the filter taps are updated without applying any transformation to the input. Once the transformation is obtained, it is applied to the rest of the data stream (TDLMS). It can be verified from Fig.3 that with DCT, the tap weight takes longer time to converge compared with PrecoG. The learning rate for LMS is static for both of these processes.

A relevant question is: how long would the initial data window be in order to compute an effective preconditioner? Table 1 in Fig. 5 and Table 2 in Fig. 6 will provide answer to that question. In each table, we have provided results on a set of attributes. First, DCT and PrecoG are tested how they perform on the estimated time-averaged R ($R_{est}$) and asymptotic R (see eq. 1). Let us take the first row of Table 1. With the initial data window of length $25$, we have obtained a transformation $U$. The DCT gives the condition number of $1.2390$ on $R_{est}$. It might appear that PrecoG performed well over the DCT if the initial $25$ data samples are considered. However, the third column reveals that with $U$ (PrecoG) at hand, it yields the condition number of $1.3550$ when the asymptotic R is considered and it is inferior to the DCT ($1.15$). This due to the fact that the $R_{est}$ is not a robust approximation of the asymptotic $R$. It is not until the data window of length $75$, we can observe that PrecoG ($1.1251$) starts performing better than DCT ($1.15$) on the asymptotic R.

Table 1 and 2 are the results on AR(1) process at two different signal correlation factors $0.95$ and $0.2$ respectively. One can see that PrecoG starts performing better than DCT with the data window length $25$ in case of $\rho=0.2$. This observation is harmonious with the theory of TDLMS. AR(1) with $\rho=0.95$ is a highly correlated datastream. Therefore, a longer data window is required to estimate $R$ in order to decorrelate the data. Fig. 4(b) shows the transformation matrices at three different initial length of data window for AR(1) process with $\rho=0.2$. When compared to Fig. 3(a) (AR(1) with $\rho=0.95$), it is evident that with smaller size data window, PrecoG contains a large number of inhibitory drives when the signal correlation factor is relatively low.

VI-B $2^{nd}$* order autoregressive process*

In Fig. 8(a), we present the condition ratios computed by applying the algorithms on the autocorrelation matrices of eight $2^{nd}$ order autoregressive process with parameters $(\rho_{1},\rho_{2})$ [10].

The input autocorrelation matrix $R_{N}$ of such process is given by $R_{N}~{}=~{}c_{1}R_{N}(\rho_{1})+c_{2}R_{N}(\rho_{2})$. $R_{N}(\rho_{1})$ and $R_{N}(\rho_{2})$ are two Toeplitz matrices, similar to $R_{N}$ of $1^{st}$ order Markov process. $c_{1}$ and $c_{2}$ are constants and are given by $c_{1}=\frac{\rho_{1}(1-\rho_{2}^{2})}{(\rho_{1}-\rho_{2})(1+\rho_{1}\rho_{2})}$, $c_{2}=\frac{-\rho_{2}(1-\rho_{1}^{2})}{(\rho_{1}-\rho_{2})(1+\rho_{1}\rho_{2})}$. As shown in Fig.8(a), PrecoG outperforms DCT in all the above cases, implying that it has better decorrelating ability than the DCT.

Next, we apply PrecoG to a simulated $3200$-length data sample from an AR(2) process with $\rho_{1}=0.6$ and $\rho_{2}=0.9$.The AR(2) data power is set to unity. We assume that there are three filter taps. The convolution filter $h=[1~{}0.8~{}-3]$ to generate the desired output. The convolved response is additively corrupted with white Gaussian noise with normalized signal power of $0.01$. For each data $x$, we set LMS learning step as $0.001$, and the power normalization factor $beta$ as $0.85$. The length of initial data window to estimate the autocorrelation matrix is set as $200$. Fig. 7(a) shows the behavior of a filter tap in the cases of DCT-LMS and PrecoG-LMS where each datasample is processed. The tap weight attains convergence faster in PrecoG accompanied by an accelerated reduction in error (Fig. 7(b)) when compared with the DCT.

VI-C Hebbian learning and Hebb-LMS algorithm

In 1949, Donald Hebb [42] postulated that concurrent synaptic changes occur as a function of pre and post synaptic activity, which was elegantly stated as ”neurons that are wired together fire together.” This classical remark of Hebb attempted to explain the plausible role of synaptic changes in learning and memory. Hebbian and LMS were predominantly regarded as two distinct forms of learning. Hebbian form of learning is unsupervised in nature, whereas LMS is primarily supervised. However, Hebbian rule, when translated into computational filters, leads to instability. Neuroscience researchers explained that neurons and the network collectively maintain stability by scaling the data at the input synapses (homeostatic plasticity) [43, 44]. From the computational point of view, a parallel of this Hebbian learning and synaptic scaling with our work is given in Appendix (Hebb-LMS algorithm). In 2019, Widrow proposed unsupervised Hebb-LMS algorithm that contains an analogue of homeostatic plasticity while using a sigmoid. PrecoG is found to show its efficacy there.

The relevance of this Hebbian learning and homeostatic plasticity (synaptic scaling) in our work can be speculated from a computational point of view. In supervised settings, the error at the output corresponding to an input data is backpropagated and eventually updates the weights. A data vector is an array of real numbers. The positive entries in a vector collectively constitute an excitatory drive and the negative entries constitute an inhibitory drive to the postsynaptic neuron, where the drives are summed. An ideal transformation (TDLMS) should have the property that it can regulate the magnitude of both the drives. A suitable transformation is expected to assess the balance of excitatory and inhibitory drives and modulate the input to achieve decorrelation. We resort to the input autocorrelation matrix that may provide such assessment. It is to remind the reader that input data is transformed prior to passing to the filter. Therefore, the transformation ‘scales’ the data for the LMS filter. So, it acts as a regularizer for the adaptation of filter weights. If this philosophy is true, PrecoG should also work in the unsupervised settings. In 2019, Widrow proposed unsupervised Hebb-LMS algorithm that contains an analogue of homeostatic plasticity while using a sigmoid. PrecoG is found to show its efficacy there.

The experimental setting considers a $3200$ length data generated by an AR(1) process with a Gaussian noise with zero mean and a variance of $3$. We assume that there are four filter taps. For each data $x$, the desired output is $Sigmoid(W^{T}x)$. We set the $slpha$ of the sigmoid as $0.5$, $gamma$ as $0.5$, the LMS learning rate as $0.001$, and power normalization factor $beta$ as $0.85$. The length of initial data window to estimate the autocorrelation matrix is set as $100$. We run $500$ epochs to inspect the behavior of error. In each epoch, weights are continuously updated after each data sample is passed and the datastream is randomly shuffled. We take the Euclidean norm of the $3200$ errors in each epoch to plot. We keep this setting fixed for all Hebb-LMS, DCT-Hebb-LMS and PrecoG-Hebb-LMS on the data. Similar conclusion can be drawn when PrecoG and DCT are tested on an AR(2) process with a Gaussian noise with zero mean and variance $5$ (fig. 9)

VII Discussion

LMS filters play pivotal roles in several applications, and convergence of the algorithm in terms of the filter weights, which are updated at each iteration, poses critical challenges to the quality of performance of LMS filters. TDLMS offers a solution to overcome the convergence issue by transforming the data prior to channeling it to the LMS filter. Conventional transforms serve as a palette of such transformation. However, the convergence of LMS is mediated by the shape of the real-life data manifold. If the shape of the data manifold is known a priori, only then we can robustly evaluate the efficacy of the transformations. This issue was largely overlooked. PrecoG provides answer to that problem by presenting an optimization framework that intelligently encodes data manifold into the transformation matrix.

Apart from symmetry by definition in case of real valued data, autocorrelation matrices containing distinct structures, such as diagonally dominant, Toeplitz, Hankel and circulant appear in special circumstances. For example, the autocorrelation matrix of $1^{st}$ order Markov process possesses Toeplitz structure. If the autocorrelation at two time points depends only on the time difference, then the autocorrelation matrix becomes Toeplitz. However, these assumptions are scarcely valid in case of data with dynamic changes over time. PrecoG evaluates the data manifold after each interval and accordingly finds the unitary preconditioner. Users can set a time interval after which the preconditioning matrix will be computed periodically to incorporate dynamic changes into the matrix. Nevertheless, this flexibility comes at the price of minor additional computational overhead. To assess the computation time, we need to look at eq.14 and Appendix.

PrecoG is shown to significantly outperform the conventional transformation matrices that are considered in this paper on autoregressive datasets and different linear systems of equations. PrecoG has been shown ( in Table 1 and Table 2 ) to have improved performance when the autocorrelation matrix of the process is estimated, making PrecoG amenable to be deployed in on-line scenarios. For autoregressive datasets, we also examined the variation of the initial window length with the signal correlation factor. It was shown to be effective in unsupervised settings also, when we investigated PrecoG’s performance in the Hebbian-LMS settings.

VIII Conclusion

In this work, we present a method to obtain a unitary split preconditioner by utilizing nonconvex optimization, graph theory and first-order perturbation theory. We demonstrate the efficacy of our approach over prevalent state-of-the-art techniques on Markov datasets and linear systems of equations. We inspect PrecoG in supervised and unsupervised settings. As a future endeavor, we will attempt to exploit the signal structure and embed this structure into the optimization framework by including a set of constraints. In continuation, we will try to extend our approach to solve a sparse underdetermined linear system of equations in order to implement dictionary learning. In addition, we will inspect the behavior of TDLMS filters by assigning the characteristics of computational neurons on the LMS taps.

Appendices

**LMS and TDLMS algorithms

**Let $X$ be a real-valued data that is channeled to a filter with $N$ tap delays, $W(n)=[w_{0},w_{1},...w(N-1)]$. The data vector at time $t=N$ at the taps can be written as $x(N)=[x_{N-1},x_{N-2},...,x_{0}]$. The filter output is given by $y(n)=W(n)^{T}x(n)$. Let the desired response be $d(n)$. Then, the error $e(n)=d(n)-y(n)=d(n)-W(n)^{T}x(n)$.

The LMS algorithm is formulated as follows.

[TABLE]

Assume that W(n) is independent of X(n) (because at the convergence, the values of W(n) will be stationary but X(n) will keep changing). Taking expectation on both sides of the above expression, we get

[TABLE]

Let us consider the term E(W(n))\big{(}1-2\mu R_{xx}\big{)}. let us assume $R_{xx}=R$ and $R_{xx}=U\Lambda U^{T}$. Therefore,

[TABLE]

In order to converge $E(W(n+1))$ to a stable value $2\mu R_{dx}$ in eq. (18), \big{(}I-2\mu\Lambda\big{)}^{n} must converge to zero at $n\to\infty$. This is possible only when $(1-2\mu\lambda_{i})\to 0\forall i$. So, the convergence will be quickly achieved if $\frac{\lambda_{max}}{\lambda_{min}}\to 1$.

In TDLMS, let the transformation be $T\in\mathcal{R}^{N\times N}$ and the transformed data vector is $x^{{}^{\prime}}(N)=Tx(N)$. The following step is the normalization of power.

[TABLE]

Next, $v(n)$ is channeled to the LMS filter. the expression can be obtained simply by replacing $x(n-i)$ with $v(n-i)$ in eq. (17).

**Hebb-LMS algorithm

**In 1949, Donald Hebb [42] postulated that concurrent synaptic changes occur as a function of pre and post synaptic activity, which was elegantly stated as ”neurons that are wired together fire together”. This classical remark of Hebb attempted to explain the plausible role of synaptic changes in learning and memory. However, it was just a conjecture at that time. Twenty years later, in a landmark paper, Bliss and Lomo [45] showed the existence of long-term potentiation (LTP) with the help of applying a train of brief, high frequency stimulation, known as tetanus, to the Hippocampus. This discovery led to a quest to mine the properties of LTP in the context of excitatory synaptic transmission [46, 47].

So, Hebbian and LMS were predominantly regarded as two distinct forms of learning. Hebbian form of learning is unsupervised in nature, whereas LMS is primarily supervised. However, when Hebb’s rule was applied to a linear filter to learn patterns in data by adaptively changing the weights, it posed a problem. If the concurrent pre and post synaptic activities strengthen over time, there is nothing to scale the activity down. Researchers attempted to introduce a ‘forgetting’ term inside the Hebb’s rule. In 1982, Oja [48] introduced a modification to stabilize the computational analogue of Hebb rule. In effect, the linear filter was able to learn the principal components of the input data. However, such modifications are ad hoc and, therefore, difficult to interpret.

In neuroscience, there exists a trenchant question: how does this intricate web of neurons, that is constantly undergoing physical changes maintains stability (homeostasis)? Researchers explained that neurons and the network collectively maintain the stability (for example, by regulating the average firing rate of a set of neurons) [49, 43]. One such mechanism that gives stability is called synaptic scaling [43, 44].

Widrow [41] proposed an algorithm that integrates Hebbian and LMS learning paradigms and it contains two equilibrium states for excitatory and inhibitory drives. The mathematical expression is given by

[TABLE]

Here $SGM$ is the polar sigmoid function. This algorithm is unsupervised because the target response of x(n) is $SGM(W(n)^{T}x(n))$. It contains two stable equilibrium points $1$ (excitatory) and $-1$ (inhibitory) on the SGM curve.

**Evaluation of \big{[}\frac{\partial E(p)}{\partial U_{N}}\big{]}

** Let, $M(\epsilon)=||U^{T}RU-(1+\epsilon)U^{T}RU\circ\mathcal{I}||_{F}^{2}$. Then,

[TABLE]

Eq. (22) on expansion gives

[TABLE]

Eq. (23) is obtained using the fact that $UU^{T}=\mathcal{I}$. By performing the partial derivative,

[TABLE]

Next, $\frac{\partial E(p)}{\partial U}$ is computed using $\frac{\partial M}{\partial U}$ as $\frac{\partial E(p)}{\partial U}=\frac{\partial M(+\epsilon_{1})}{\partial U}+\frac{\partial M(-\epsilon_{2})}{\partial U}$.

[TABLE]

[TABLE]

**Proof of eq. (15)

**Let us assume $A_{0}$ is a real-valued, finite-dimensional ($\in\mathcal{R}^{N}$), positive semi-definite matrix with eigenvectors $U_{0}=[\boldsymbol{u}_{1},\boldsymbol{u}_{2},...,\boldsymbol{u}_{N}]$ and the corresponding eigenvalues $\lambda_{0}=[\lambda_{1},\lambda_{2},...\lambda_{N}]$ (arranged in the decreasing order). In our case, $A_{0}=L$, the graph Laplacian matrix. By definition $\lambda_{N}=0$.

Let us also assume that $A$, $U$ and $\lambda$ are real-valued functions of a continuous parameter $\tau$ and they are analytic in the neighborhood of $\tau_{0}$ such that $A(\tau_{0})=A_{0}$, $U(\tau_{0})=U_{0}$ and $\lambda(\tau_{0})=\lambda_{0}$.

[TABLE]

Now, because of the facts that $\boldsymbol{u}_{i}$ is analytic and $||\boldsymbol{u}_{i}||=1$ (orthonormal), $\dot{\boldsymbol{u}}_{i}\bot\boldsymbol{u}_{i}$. Taking inner product with $\boldsymbol{u}_{i}$ on both sides and plugging \big{<}\dot{\boldsymbol{u}}_{i},\boldsymbol{u}_{i}\big{>}=0, it can be found that

[TABLE]

Here, \big{<}A\dot{\boldsymbol{u}}_{i},\boldsymbol{u}_{i}\big{>}=\big{<}\dot{\boldsymbol{u}}_{i},A\boldsymbol{u}_{i}\big{>} (because $A$ is symmetric (self-adjoint)). Using $A\boldsymbol{u}_{i}=\lambda_{i}\boldsymbol{u}_{i}$, we get \big{<}\dot{\boldsymbol{u}}_{i},A\boldsymbol{u}_{i}\big{>} = \big{<}\dot{\boldsymbol{u}}_{i},\lambda_{i}\boldsymbol{u}_{i}\big{>}= \lambda_{i}\big{<}\dot{\boldsymbol{u}}_{i},\boldsymbol{u}_{i}\big{>}=0. From eq. Appendices, we obtain the expression for $\dot{\lambda}_{i}$ as \dot{\lambda}_{i}=\big{<}\dot{A}\boldsymbol{u}_{i},\boldsymbol{u}_{i}\big{>}.

Taking the inner product of eq. Appendices with $\boldsymbol{u}_{j};~{}j\neq i$ and using the fact that $\boldsymbol{u}_{i}\bot\boldsymbol{u}_{j}$ we obtain

[TABLE]

Let us consider $\tau=a_{mn}$, where $a_{mn}$ is an entry of $A$.

[TABLE]

It can be seen that $\dot{\boldsymbol{u}}_{i}$ is spanned by the eigenvectors of $A$, when all the eigenvalues are numerically different. Due to the fact that $A$ is symmetric, $\frac{\partial A}{\partial a_{mn}}$ is also symmetric. Therefore, eq. 28 can be rearranged as

[TABLE]

The above expression is valid for $\lambda_{i}\neq\lambda_{j}$. We can use the fourth step of eq. Appendices to derive the expression of $\dot{\boldsymbol{u}}_{i}$ in case of $\lambda_{i}=\lambda_{j}$ (multiplicity of eigenvalues).

[TABLE]

Combining eq. Appendices and Appendices,

[TABLE]

Here $\delta$ is the Kronecker delta function.

**Proof of eq.(16)

**By definition,

[TABLE]

where $L_{mn}$ indicates the $(m,n)^{th}$ element of the symmetric matrix $L$. Let us first find the expression for $\lambda_{i}\neq\lambda_{j}$. Inserting eq. Appendices by replacing $A$ with $L$ to the above equation, we get

[TABLE]

From eq. 33, it can be observed that the eigenvector $\boldsymbol{u}_{q}$ is projected by the operator $\frac{1}{(\lambda_{k}-\lambda_{q})}B\frac{\partial W}{\partial w_{j}}B^{T}$ prior to the inner product.

Applying the same definition the final expression can be found by using eq. 31.

[TABLE]

Proof: $\frac{\partial L}{\partial u_{kl}}=U_{N}\Gamma J^{mn}+J^{nm}\Gamma U_{N}^{T}$ is non-invertible

By definition, $J^{mn}=\delta_{mk}\delta_{np}$. Let $F=\Gamma J^{mn}$. Then

[TABLE]

Therefore,

[TABLE]

So, $U_{N}\Gamma J^{mn}$ has one column (at $y=p$) that has a set of non-zero entries. It is evident that $H=(U_{N}\Gamma J^{mn})^{T}=J^{nm}\Gamma U_{N}^{T}$. Therefore, $\frac{\partial L}{\partial u_{kl}}$ has one row and one column of possibly non-zero entries by construction. Assume $m$ is that row index and the corresponding entry of $H$ is $[H_{m1},H_{m2},...,H_{mN}]$. Assume that $q\in\{1,2,...,N\}$ is the column index, where $H$ receives entries from $U_{N}\Gamma J^{mn}$. One can find any two columns $i,j\in\{1,2,...,N\}-\{q\}$ such that the following transformation $C_{i}\longrightarrow\frac{C_{i}}{H_{mi}}$ and $C_{j}\longrightarrow\frac{C_{j}}{H_{mj}}$ would give two identical columns.

**Complexity

**

The update of weights $w_{i}$ requires the computation of three partial derivatives (see Appendix). *Using eq. 34

*Note that $B\frac{\partial W}{\partial w_{j}}B^{T}$ is fixed for $w_{j}$. By definition, $L=BWB^{T}$. So, $E=\frac{\partial L}{\partial w_{j}}=B\frac{\partial W}{\partial w_{j}}B^{T}$. For a simple, connected graph, there is no self-loop. Let us assume $w_{j}=w_{kq}$, meaning that the edge between the $k^{th}$ and $q^{th}$ vertices has a weight $w_{j}$. $L$ is a linear function of $w_{j}$. Then, $E$ has exactly four non-zero entries - $E(k,k)=E(q,q)=1$, $E(k,q)=E(q,k)=-1$. Therefore, $B\frac{\partial W}{\partial w_{j}}B^{T}\boldsymbol{u}_{p}$ requires constant computation. For example, if $w_{j}=w_{13}$, then $B\frac{\partial W}{\partial w_{j}}B^{T}\boldsymbol{u}_{p}$= $[u_{p1}-u_{p3},0,u_{p3}-u_{p1},0,0,..]^{T}$, having non-zer entries at the $1^{st}$ and $3^{rd}$ locations. It implies that \big{<}\boldsymbol{u}_{i},B\frac{\partial W}{\partial w_{j}}B^{T}\boldsymbol{u}_{q}\big{>} requires constant computation. Now, each $\boldsymbol{u}_{i}$ has length $N$. For each $i$, $\frac{\partial\boldsymbol{u}_{i}}{\partial w_{j}}$ requires $\mathcal{O}(N)$ time to add for all $(N-1)$ vectors. So, the time complexity to compute $\frac{\partial U_{N}}{\partial w_{j}}$ is $\mathcal{O}(N^{2})$, which is a significant improvement over the last scheme. The improvement is due to the non-existence of matrix inversion.

It is to note that the constructions of preconditioners for solving linear systems by comparative methods such as, Jacobi ($\mathcal{O}(N^{2})$), successive over-relaxation (SOR) ($\mathcal{O}(N^{3})$), symmetric SOR ($\mathcal{O}(N^{3})$), Gauss-Seidel ($\mathcal{O}(N^{3})$) have faster associated run times compared to PrecoG. Here, complexity accounts for the inversion of each preconditioner matrix. However, the acceleration in the convergence of the LMS filter using PrecoG is also expected to partially compensate for the computational overload of PrecoG.