Tensor network square root Kalman filter for online Gaussian process regression
Clara Menzen Manon Kok Kim Batselier



 
Delft]Clara Menzencm.menzen@gmail.com,
Delft]Manon Kokm.kok-1@tudelft.nl,
Delft]Kim Batselierk.batselier@tudelft.nl

[Delft]Delft Center for Systems and Control, Delft University of Technology, Mekelweg 2, 2628 CD, Delft, the Netherlands

          
                           
Square root Kalman filtering, tensor network, Gaussian processes, recursive estimation.
                             
                                        
§ ABSTRACT 
The state-of-the-art tensor network Kalman filter lifts the curse of dimensionality for high-dimensional recursive estimation problems.
However, the required rounding operation can cause filter divergence due to the loss of positive definiteness of covariance matrices.
We solve this issue by developing, for the first time, a tensor network square root Kalman filter, and apply it to high-dimensional online Gaussian process regression.
In our experiments, we demonstrate that our method is equivalent to the conventional Kalman filter when choosing a full-rank tensor network.
Furthermore, we apply our method to a real-life system identification problem where we estimate 4^14 parameters on a standard laptop.
The estimated model outperforms the state-of-the-art tensor network Kalman filter in terms of prediction accuracy and uncertainty quantification.





§ INTRODUCTION

In a time when data-driven AI models are trained on an exponentially growing amount of data, it is crucial that the models can be adapted to newly observed data without retraining from scratch.
These online or recursive settings are present in many fields including system identification <cit.>, sensor fusion <cit.>, robotics <cit.>, and machine learning <cit.>.

While Bayesian algorithms, like widely-used Gaussian processes (GPs)  <cit.> are well-suited for an online setting, they are associated with potentially high computational costs.
Standard GP regression using a batch of N observations has a cubic cost in N, i.e., 𝒪(N^3).
The number of observations is growing in an online setting, so the cost increases each time step and can become a computational bottleneck.

There are numerous parametric approximations to address scalability in batch settings, including sparse GPs <cit.> and reduced-rank GPs <cit.>, which both have a complexity of 𝒪(NM^2), M being the number of inducing inputs and basis function for the respective method.
Structured kernel interpolation for sparse GPs <cit.> reduces the complexity further to 𝒪(N+DM^1+1/D), D being the number of input dimensions.

Parametric approximations allow for a straightforward recursive update, where the posterior distribution from the previous time step is used as a prior for the current time step <cit.>.
In this context, online GPs have been used, e.g.,  for GP state-space models <cit.>, rank-reduced Kalman filtering <cit.> and recursive sparse GPs <cit.>.

In this paper, we consider the online parametric GP model given by

    y_t =  ϕ   (_t)^⊤_t + ϵ_t,  ϵ_t∼𝒩(0,σ_y^2),
       _t-1∼𝒩(_t-1,𝐏_t-1),

where y_t is a scalar observation at discrete time t, ϕ(·) are basis functions that map a D dimensional input vector _t to a feature space, _t∈ℝ^M are the parameters at time t, and σ_y^2 denotes the variance of the measurement noise ϵ_t which is assumed to be i.d.d. and zero-mean Gaussian.
With (<ref>), the posterior distribution  is computed each time step using the estimate _t-1 and covariance matrix _t-1 of the previous time step as a prior. 

We consider commonly used product kernels with a feature map given by

    ϕ(_t) = ϕ^(1)(_t) ⊗⋯⊗ϕ^(d)(_t) ⊗⋯⊗ϕ^(D)(_t),

where ϕ^(d)(_t)∈ℝ^I with I being the number of basis functions in the dth dimension.
The resulting number of basis functions is M=I^D, growing exponentially with the input dimension D.
Several tensor network (TN)- based methods have been proposed to break this curse of dimensionality and achieve a linear computational complexity in D.
In the batch setting, <cit.> and <cit.> give solutions for the squared exponential and polynomial kernel, respectively.
In the online setting, the state-of-the-art method is the tensor network Kalman filter (TNKF) <cit.>, where the Kalman filter time and measurement update are implemented in TN format. 

While the TNKF lifts the curse of dimensionality, it has a significant drawback.
The TNKF requires a TN-specific rounding operation <cit.>, which can result in covariance update losing positive   definiteness <cit.>, resulting in the divergence of the filter.



This paper resolves this issue by computing the square root covariance factor in tensor train (TT) format instead.
Our approximation represents the M× M square root covariance factor as a tensor train matrix (TTm). 
This is motivated by prior square root covariance factors of product kernels having a Kronecker product structure, which corresponds to a rank-1 TTm.
In addition, work by <cit.> and <cit.> approximate the covariance matrix as a rank-1 TTm.
This work generalizes the rank-1 approximation to higher ranks which results in better prediction accuracy and uncertainty quantification.
We call our method the tensor network square root Kalman filter (TNSRKF).

We show in experiments that the TNSRKF is equivalent to the standard Kalman filter when choosing full-rank TTs.
In addition, we show how different choices of TT-ranks affect the performance of our method.
Finally, we compare the TNSRKF to the TNKF in a real life system identification problem with 4^14 parameters, and show that contrary to the TNKF, our method does not diverge.




§ PROBLEM FORMULATION


Similar to the TNKF, we build on standard equations for the measurement update of the Kalman filter given by

    _t    = ϕ_t^⊤_t-1ϕ_t + σ_y^2
    _t    = _t-1ϕ_t _t^-1
    _t    = _t-1 + _t(y_t - ϕ_t^⊤_t-1)
    _t    = (𝕀_M-_tϕ_t^⊤)_t-1(𝕀_M-_tϕ_t^⊤)^⊤ + σ_y^2_t^⊤_t,

where _t denotes the innovation covariance and _t denotes the Kalman gain. 
Note that for a scalar measurement, _t is a scalar and _t a vector, whereas in the case of multiple measurements per time step, they are matrices.
In this way, we recursively update the posterior distribution of the parametric weights from (<ref>), i.e.,  p(_t|_1:t).
For product kernels with a feature map given in (<ref>), it is _t∈ℝ^I^D and _t∈ℝ^I^D× I^D.
In this case, the Kalman filter suffers from the curse of dimensionality.

The first tensor-based Kalman filter, the TNKF <cit.>, solved the curse of dimensionality and implements (<ref>)-(<ref>) in TT format, where the weights are represented as a TT and the covariance matrix as a TTm.
During the updates, the algebraic operations in TT format increase the TT-ranks of the involved variables, according to <cit.>.
To counteract the rank increase and keep the algorithm efficient, the TNKF requires an additional step called TT-rounding <cit.>.
This SVD-based operation transforms the TT or TTm to ones with smaller TT-ranks.
TT-rounding can result, however, in the loss of positive (semi-) definiteness.

To avoid this issue, we implement the square root formulation of the Kalman filter (SRKF), as described e.g. in <cit.>, in TT format.
The SRKF expresses (<ref>)-(<ref>) in terms of the square root covariance factor, given by

    _t = [ (𝕀_M-_tϕ_t^⊤)_t-1             σ_y_t ],





which consists of a concatenating matrices that increases the number of columns.
For the next update, _t needs to be transformed back to its original size.
In the SRKF, this is done by computing a thin QR-decomposition <cit.> of _t given by

    _t^⊤_(M+1)× M = 𝐐_t_(M+1)× M𝐑_t_M× M

and replacing _t by 𝐑_t.
The orthogonal 𝐐_t-factor can be discarded since

    _t = _t_t^⊤ = 𝐑_t^⊤𝐐_t^⊤𝐐_t_𝐈_M𝐑_t = 𝐑_t^⊤𝐑_t.


In TT format, performing the QR-decomposition as in (<ref>) is not possible.
We solve this issue by proposing an SVD-based algorithm in TT format that truncates _t back to its original size.
















§ BACKGROUND ON TENSOR NETWORKS



 §.§ Tensor networks

Tensor networks (TNs), also called tensor decompositions, are an extension of matrix decompositions to higher dimensions.

In this paper, we use a specific architecture of TNs, called TTs <cit.> to approximate the weight vector's mean as discussed in Section <ref>, and a TT matrix (TTm) <cit.> to approximate the square root covariance factor, as discussed in Section <ref>.

In this context, we denote TTs representing vectors as a lower-case bold letter, e.g. _t, and their components, called TT-cores, as capital calligraphic bold letters, e.g. W^(d).
TT matrices are denoted by upper-case bold letters, e.g. _t and their corresponding TTm-cores as capital calligraphic bold letters, e.g. L^(d).





  §.§.§ Tensor train vectors


As depicted in Fig. <ref>(a), a TT vector consists of interconnected three-way tensors, called TT-cores, visualized as nodes with three edges.
Each edge corresponds to an index of a TT-core and connected edges are summations over the involved indices.
Each TT-core is connected by two edges, called TT-ranks, to their neighbouring TT-cores, except for the first and last TT-core, which outer TT-ranks are by definition equal to one.

For the purpose of this paper, consider a TT that represents the mean of the weight vector _t∈ℝ^M.
The TT-cores, denoted by W_t^(1),⋯,W_t^(d),⋯,W_t^(D) are each of size R_d× I× R_d+1, where R_d and R_d+1 are the TT-ranks and I is the size of the non-connected edge such that .
By definition R_1=R_D+1=1.
Without the loss of generality, we use TT-cores with equal TT-ranks R_.
The storage complexity of _t without TNs is 𝒪(I^D) and in TT format 𝒪(DIR_^2), where lower TT-ranks R_ will result in more efficient representations.

An important characteristic of a TT for numerical stability is that it can be transformed into the site-d-mixed canonical format.
Site-d-mixed canonical format <cit.>
    A TT _t in site-d-mixed canonical format is given by 
    
    _t=𝐆_d,t_t^(d),

    where 𝐆_d,t∈ℝ^M× R_ IR_ is an orthogonal matrix computed from all TT-cores except the dth and  is the vectorization of the dth TT-core.
    In this format, the TT representation is linear in the dth TT-core when all other TT-cores are fixed. 




  §.§.§ TT matrices and tall TT matrices


A TTm consists of interconnected four-way tensors, as depicted in Fig. <ref>(b).
Analogous to the TT, the TTm components and connected edges are called TTm-cores and TTm-ranks, respectively, where each TTm-core has two free edges, the row and column indices.

For the purpose of this paper, consider a TTm representation of the square root covariance factor _t∈ℝ^M× M.
The TTm-cores L_t^(1),⋯,L_t^(d),⋯,L_t^(D) are of size R_d× I× J × R_d+1, where I and J are the row and column indices, indicated in Fig. <ref>(b) as red and and blue edges respectively, such that  and M=J^D.
By definition, R_1=R_D+1=1, and for this paper, we generally assume that all other TTm-ranks R_2=⋯=R_D=R_ are equal.
The storage complexity of _t without TNs is 𝒪(I^D× I^D) and in TTm format 𝒪(DR_^2 IJ).

A TTm can also be written in terms of the site-d-mixed canonical format as defined in Definition <ref>, but it requires to be transformed into a TT first.
This can be done by combining the row and column indexes into one index, which represents a kind of vectorization of the matrix represented by the TTm. 
Note, however, that the indices are not ordered as in conventional vectorization.
A site-d-mixed canonical format of a TTm is given by 

    vec(_t) =𝐇_d,t𝐥_t^(d),

where 𝐇_d,t∈ℝ^2M× R_ IJR_ is computed from all the TTm-cores but the dth, and 𝐥^(d)_t∈ℝ^R_ IJR_.

To recompute _t back in its original like in the QR step of the SRKF (see (<ref>)), we need a special of a TTm, the tall TTm, as well as a thin SVD in TTm format.
Tall TTm <cit.>
    A tall TTm, as depicted in Figure <ref>(c), has only one TTm-core with both a row and column index, while all other TTm-cores have only row indices.
    Then, the TTm represents a tall matrix with many more rows than columns.


Thin SVD in TTm format <cit.>

   Consider a TTm in site-d-mixed canonical format, where the dth TTm-core is the one that has the column index, .
   The SVD of L^(d) reshaped and permuted in to a matrix of size R_ IR_× J, is given by
   
    𝐔^(d)𝐒^(d)(𝐕^(d))^⊤.

    Now replace the dth TTm-core by 𝐔^(d) reshaped and permuted back to the original TTm-core dimensions.  
    
   Then the thin SVD is given by the TTm with the replaced TT-core as the orthogonal 𝐔-factor, and 𝐒^(d)(𝐕^(d))^⊤ as the 𝐒𝐕^⊤-factors, as depicted in Fig. <ref>(d).





§ TENSOR-NETWORKED SRKF

We propose our method combining efficient TN methods with the SRKF formulation for online GP regression.
More specifically, we recursively compute the posterior distribution of the parametric weights in (<ref>) from the measurement update of the Kalman filter.
To achieve this, we update the mean _t∈ℝ^M as a TT (Section <ref>), and the square root factor _t∈ℝ^M× M as a TTm (Section <ref>).

All computations are summarized in Algorithm <ref>, which outputs the posterior weight distributions p(_t|_1:t)=𝒩(_t,_t), and the predictive distributions for a test input p(f_*,t)=𝒩(m_*,t,σ_*,t^2).



 §.§ Update of weight mean


For updating the mean with a new measurement y_t∈ℝ, we compute (<ref>) in TT format.
In the original tensor-based KF <cit.>, the two terms in equation (<ref>) are summed together in TT format, which increases the TT-ranks.
To avoid this rank increase, we apply a commonly-used optimization algorithm from the tensor community, called the alternating linear scheme (ALS) <cit.>.
The ALS computes a TT by updating one TT-core at a time while keeping all other TT-cores fixed.
The optimization problem to be solved is given by

    _tmin   _t-1 + _t(y_t-ϕ_t^⊤_t-1) - _t _F^2
       s.t.  _t  being a low-rank TT,

where the subscript F stands for the Frobenius norm, and _t-1 is the estimate from the last time step playing now the role of the prior for the current time step.

Inserting (<ref>) in (<ref>), thus making use of the site-d-mixed canonical format from Definition <ref>, gives the optimization problem for the update of one TT-core 

    _t^(d)min𝐆_d,t^⊤(_t-1 + 𝐊_t(y_t-ϕ_t^⊤_t-1)) - _t^(d)_F^2.

In one so-called sweep of the ALS, (<ref>) is solved for each TT-core once.
A stopping criterion for the convergence of the residual in (<ref>) determines the total number of sweeps.




 §.§ Update of square root covariance factor


To compute the covariance matrix with the standard covariance update in the measurement update, see (<ref>), we recursively compute the square root covariance factor _t as defined in (<ref>) such that _t=_t_t^⊤.
To achieve this, we use the ALS to solve (<ref>) (ALS step) and then we transform _t as in (<ref>) back to its original size (QR step).



  
*ALS step
In this step, we use the ALS to compute a TTm representing _t. 
We solve the optimization problem given by

    min__t   [ (𝕀_M-_tϕ_t^⊤)_t-1             σ_y_t ] - _t _F^2 
       s.t.  _t  being a low-rank TTm,

where _t-1 is the estimated square root covariance factor from time step t-1 now serving as the prior.
The original ALS algorithm is defined for TTs, so we must adapt it for TT matrices.

For this, it is necessary to use the site-d-mixed canonical form for TT matrices, as described in Section <ref> above (<ref>).
In addition, we need to horizontally concatenate two matrices in TTm format, which can be done by summing two matrices of size M× 2M such that (<ref>) becomes

    min_𝐥_t^(d)𝐇_d,t^⊤vec([ 1 0 ]⊗ (𝕀_M-_tϕ_t^⊤)_t-1) . 
    .+𝐇_d,t^⊤vec([ 0 1 ]⊗[     1 0_M-1 ]⊗σ_y_t)   - 𝐥_t^(d)) _F^2,

where vec denotes the vectorizations of the involved TT matrices.



  
*QR step
The optimization problem given by (<ref>) requires concatenating a matrix with a column vector.
In TT format, this results in a TTm of size M× 2M.
For the TTm-cores of _t this means that one TTm-core, we call it the augmented core, is of size R_× I× 2J× R_.
Before serving as a prior for the next time step, a QR step as in (<ref>) is required to transform _t back to its original size.
We use an SVD-based algorithm in TN format to transform _t of size M× 2M back to M× M, as described in Algorithm <ref>.



 §.§ Predictions

To perform GP predictions we compute the predictive distribution for a test output f_*,t=ϕ(_*)^⊤_t with mean and variance given by

    m_*,t   = ϕ(_*)^⊤_t
    σ^2_*,t   = ϕ(_*)^⊤_t_t^⊤ϕ(_*).


Given _t as a TT and _t as a TTm, we can compute (<ref>) directly in TN format without explicitly reconstructing the mean vector and square root factor.
For a test input _*, Fig. <ref> illustrates the computation of (a) the predictive mean m_*,t, (b) the predictive covariance σ^2_*,t.















§ IMPLEMENTATION


In this section, we give a detailed description of the non-straight-forward TN operations to update the mean estimate _t and square root covriance factor _t as described in Algorithm <ref>.
The leading complexities of the mean and square root covriance factor update are given in Table <ref>.








 §.§ Updating _t in TN format

In the following sections, we discuss the implementation of (<ref>) for the mean update (Algorithm <ref>, line <ref>), and we describe how the mean is initialized in TT format (Algorithm <ref>, line <ref>).



  §.§.§ Implementation of 𝐆_d,t^⊤(_t-1 + 𝐊_t(y_t-ϕ_t^⊤_t-1))

To compute the TT representing the mean estimate _t, we implement the ALS to solve (<ref>) (

The following example illustrates the update of one TT-core during the ALS.
TT-core update with ALS
    Take a D=5 dimensional weight vector in TT format with M_1=M_2=M_3=M_4=M_5=10 basis functions in each dimension, resulting in 10^5 parameters and uniform TT-ranks of R_2=R_3=R_4=4.
    Say, we are currently updating the third TT-core W_t^(3)∈ℝ^4×10× 4 using
    
    _t^(3)_160×1 = 𝐆_3,t^⊤_160× 10^5(_t-1_10^5× 1 + 𝐊_t_10^5× 1(y_t-ϕ_t^⊤_t-1)_1×1).

    We first multiply over the large dimension of 10^5 in 𝐆_3,t^⊤_t-1 and 𝐆_3,t^⊤_t(y_t-ϕ_t^⊤_t-1). 
    In TT format, this matrix-vector-multiplication is done core by core, thus avoiding 
    the explicit multiplication. 
    Finally, we sum two vectors of size 160.
    

Figure <ref> illustrates the multiplication of 𝐆_d,t^⊤_t-1 in TT format, resulting in a tensor of size R_× I× R_.
The multiplication of between 𝐆_d,t^⊤ and _t(y_t-ϕ_t^⊤_t-1) works in the same way after firstly computing ϕ_t^⊤_t-1 in TN format and secondly multiplying one arbitrary TT-core of _t by the scalar (y_t-ϕ_t^⊤_t-1).

During the update of the dth TT-core, the TT is in site-d-mixed canonical format.
Before updating the next TT-core, either the (d-1)th or the (d+1)th, the site-(d-1)-mixed or site-(d+1)-mixed canonical format is computed.
Note that because of the recursive property, updating every TT-core once with a new measurement is usually sufficient for the residual of (<ref>) to converge.






  §.§.§ Initialization of _0 and _1

For the first time step t=1 of Algorithm <ref>, we choose a zero-mean assumption for the prior estimate _0.
The following Lemma explains how this can be implemented in TT format.

Zero-mean prior in TT format <cit.>
    Consider a vector with all entries equal to zero.
    In TT format, such a vector is given by a TT in site-d-mixed canonical format, where the dth TT-core contains only zeros. 


In addition, Algorithm <ref> requires an initial guess for _1 to compute 𝐆_d,1 from all TT-cores of _1, except the dth.
For this, we set _1=_0.





 §.§ Updating _t in TT format



To compute the TTm representing _t, we implement the ALS to solve (<ref>) (Algorithm <ref>, line <ref>).
The following example illustrates the update of one TTm-core during the ALS.

TTm-core update with ALS
    Take a D=5 dimensional TTm representing _t∈ℝ^M× M, where we are currently updating the third TTm-core. 
    We have I=10 and J=10, where the third TTm-core is augmented, and R_=4.
    We update L_t^(3)∈ℝ^4×10×20× 4 using
    
    𝐥_t^(3)_3200×1 =   𝐇_d,t^⊤_3200×2·10^10   vec([ 1 0 ]_1× 2⊗_t-1_10^5× 10^5)
    
            - 𝐇_d,t^⊤_3200×2·10^10   vec([ 1 0 ]_1× 2⊗_tϕ_t^⊤_t-1_10^5× 10^5)
    
            + 𝐇_d,t^⊤_3200×2·10^10   vec([ 0 1 ]_1× 2⊗[     1 0_M-1 ]_1× 10^5⊗σ_y_t_10^5× 1).

    We first multiply over the large dimension of 2·10^10 in TT format, then sum the three terms of size 3200× 1.

From Example <ref>, it follows that the three terms of (<ref>) need to be implemented.
We discuss them separately in the following sections.
We distinguish between the update of the augmented TTm-core from all other ones, which result in TTm-cores of size R_× I× 2J × R_ and R_× I× J× R_, respectively.
In the tensor diagrams (Fig. <ref>-<ref>), we depict the update for the augmented TTm-core.

Before diving in, recall from (<ref>) that 𝐇_d,t is computed from TTm-cores of _t, except the dth, where row and column indices are combined.
In the tensor diagrams, the indices are depicted not as combined because, in practice, they are generally summed over separately.
However, the vectorized format is necessary for writing down the equations in matrix form.



  §.§.§ Implementation of first term of (<ref>)

Fig. <ref> illustrates the computation of the augmented TTm-core in the first term of (<ref>), given by

    𝐇_d,t^⊤vec([ 1 0 ]⊗σ_y_t-1).

The column indices of _t-1 are indicated by the round edges that are connected to the row indices of 𝐇_d,t^⊤.
The edge containing 𝐞_1=[1  0] is connected to the dth TTm core of _t with a rank-1 connection, which corresponds to the Kronecker product in (<ref>).
The summation over the vertical and curved indices has the leading computational complexity of 𝒪(R_^4IJ) per dimension.
When updating all TTm-cores except the augmented TTm-core, the additional index of size 2 is summed over resulting in a tensor of size R_× I× J× R_.





  §.§.§ Implementation of second term of (<ref>)

Fig. <ref> illustrates the computation of 
    
    𝐇_d,t^⊤vec([ 1 0 ]⊗_t-1_t-1^⊤ϕ_t_t^-1ϕ_t^⊤_t-1),

which directly follows from the second term of (<ref>).
As shown, the row and column indices of 𝐇_d,t^⊤ are connected separately to the column and row indices of two TT matrices for _t-1, respectively.
Like in the previous term, the edge containing 𝐞_1=[1  0] is connected to the augmented TTm-core of _t with a rank-1 connection, which corresponds to the Kronecker product in (<ref>).
The leading computational complexity of 𝒪(R_^4IJ) per dimension comes from the summation over the vertical indices in the red or blue box indicated in the figure.
The most efficient order of doing the computations in Fig. <ref> was found with the visual tensor network software by <cit.>.






  §.§.§ Implementation of third term of (<ref>)

Fig. <ref> illustrates the computation of 

    𝐇_d,t^⊤vec([ 0 1 ]⊗[     1 0_M-1 ]⊗σ_y_t_t^⊤ϕ_t_t^-1),

which directly follows from the third term of (<ref>).
The row of nodes each filled with 𝐞_1=[1  0_J-1] corresponds to [1  0_M-1] from (<ref>) and their rank-1 connections to the nodes above is the second Kronecker product in (<ref>), which is done dimension-wise in TT format.
The node with 𝐞_2 corresponds to [0  1] from (<ref>) and its rank-1 connection is the first Kronecker product in (<ref>). 
The summation over the vertical indices is the leading computational complexity of 𝒪(R_^4IJ) per dimension.








  §.§.§ SVD-based QR step in TTm format


When computing (<ref>), we double the number of columns of _t compared to _t-1.
For the next time step, however, we need to transform _t back to its original size (Algorithm <ref>, line <ref>), otherwise its column size will grow with the iterations and slow down the algorithm.
The QR step, as in (<ref>), computes a full QR decomposition of _t, which cannot be done in TT format.
Instead, we compute a thin SVD in TTm format (Definition <ref>) of _t transformed into a tall TTm (here also denoted by _t) with all row indices of size IJ, except the dth which is of size I, and the dth column index of size 2J.
The J-truncated SVD of _t is then given by

    _t_MJ^D-1×2J≈𝐔_t𝐒_t_MJ^D-1× J𝐕_t^⊤_J×2J,

where 𝐔_t𝐒_t is the new _t and 𝐕_t^⊤ can be discarded because of (<ref>).
In practice, we compute an SVD of the augmented TTm-core and truncate it back to the size of R_× I× J× R_.

There is a way to make (<ref>) exact.
This is possible if the augmented TTm-core is of size R_× I× 2J R_^2× R_.
In this case, the SVD computed of the augmented TTm-core results in a square 𝐔-factor.
Since the number of columns is doubled every measurement update, the QR step can be skipped p times until 2^p=2R_^2.
Choosing smaller values for p reduces computational complexity at the cost of accuracy.

The SVD-based QR step is described in Algorithm <ref>.
The SVD of the reshaped and permuted augmented TTm-core is truncated for 2^p<2R_^2 and exact for 2^p≥ 2R_^2.





  §.§.§ Initialization of _0 and _1

At time t=1, Algorithm <ref> requires the prior square root covariance factor _0 in TTm format.
We are considering product kernels that have priors in Kronecker format.
The following Lemma describes how these types of priors can be transformed into a TTm for _0.


(Prior covariance with Kronecker structure into TTm, follows from <cit.>)
Given a prior covariance _0=_0^(1)⊗_0^(2)⊗⋯⊗_0^(D), the prior square root covariance in TTm format is given by a TTm with all ranks equal to 1, where the cores are given by _0^(1),_0^(2),…,_0^(D), each reshaped into a 4-way tensor of size 1× I× J× 1.



In addition, Algorithm <ref> requires an initial guess in TTm format for _1∈ℝ^M× 2M.
We cannot set _1 = _0 since the prior has TTm-ranks equal to one, and we may want higher TTm-ranks for _t. 
This is because the choice of the TTm-ranks of _1 determines the rank manifold on which the TTm-cores will be optimized.
We initialize the TTm-cores as random samples from a zero-mean Gaussian distribution and transform the TTm into site-d-mixed canonical format, where d is the augmented TTm-core.


§ EXPERIMENTS

In this section, we show how our method works in practice by performing online GP regression on synthetic and real-life data sets.
We evaluate our predictions based on the root mean square error (RMSE) for the accuracy of the mean and negative log-likelihood (NLL) for the uncertainty estimation.
The metrics after t measurement updates are defined as

    (RMSE)_t    = √(∑_i=1^N_*(m_*,t,i-y_*,i)^2/N_*)   and
    
            (NLL)_t    = 0.5∑_i=1^N_*log (2πσ_*,t,i^2) + (m_*,t,i-y_*,i)^2/σ_*,t,i^2,

where y_*,i is the ith measurement from the test set, m_*,t,i and σ_*,t,i are the predictive mean and variance for the ith test point, and N_* is the number of test points.

First, we show the equivalence of the full-rank TNSRKF and the conventional Kalman filter.
Then we show in a synthetic experiment how the choice of R_ and R_ impact the accuracy of the approximation.
Finally, we compare our method to the TNKF on a benchmark data set for nonlinear system identification.

All experiments were performed on an 11th Gen Intel(R) Core(TM) i7 processor running at 3.00 GHz with 16 GB RAM.
For reproducibility of the method and the experiments, the code written in Julia programming language is freely available at .



 §.§ Equivalence of full-rank TNSRKF and Kalman filter


In the first experiment, we show in which case our method is equivalent to the measurement update of the conventional Kalman filter.
We generate D=3 dimensional synthetic data sampled from a reduced-rank GP by <cit.> with a squared exponential kernel and use I=4 basis functions per dimension, such that _t∈ℝ^64×64.
The input data lies in a cuboid given by [-1  1]×[-1  1]×[-1  1] and N,N_*=100

Table <ref> shows the RMSE and NLL for test data at time t=N for different choices of p.
The TNSRKF is equivalent to the Kalman filter when both R_ and R_ are full-rank.
In addition, p must be chosen, such that the QR step, discussed in Section <ref>, is exact.
For settings with lower values for p, the method trade-in accuracy.

In the following sections, we look at scenarios where the Kalman filter can no longer be computed on a conventional laptop because both storage and computational time become unfeasible.



 §.§ Influence of the ranks on the approximation

The choice of the TT- and TTm-ranks is not obvious and can be intricate.
However, the computational budget often determines how high the ranks can be chosen.
In this experiment, we use our method to make online GP predictions on synthetic data while varying the TTm-ranks of _t, as well as the TT-ranks of _t.

We consider the Volterra kernel, a popular choice for nonlinear system identification.
It is known that the truncated Volterra series suffers from the curse of dimensionality, which was lifted in a TN setting by <cit.>.
With the notation of this paper, the basis functions ϕ of parametric model (<ref>) are a combination of monomials computed from the input sequence of the given problem.
We generate synthetic training and testing data as described in <cit.>, where D=7 and I=4 such that the number of parameters is 4^7=16384.
We set the SNR to 60, so there is relatively little noise.

Fig. <ref> shows the RMSE and NLL on the testing data for R_=2,4 and R_=2,4 over time iterations of the TNSRKF.
At t=N, the RMSE is lower for R_=2 and R_=4 than for R_=4 and R_=4.
Thus it seems that a lower value for the mean estimate represents the data better.
The NLL is the lowest for R_=4 and R_=4, which is close to the NLL for R_=4 and R_=2.
Note that the NLL for the same R_ is different for the two settings of R_, because the NLL also depends on the difference between predicted and actual measurements, thus on the accuracy of _t.




This experiment showed that the choice of R_ and R_ influence the performance of the TNSRKF.
Since higher values for the ranks also increase the computational complexity, the computational budget will determine the higher limit for the ranks.
In addition, an assumption with lower ranks may be fitting the data better in some cases.



 §.§ Comparison to TNKF for cascaded tanks benchmark data set

In this experiment, we compare our method to the TNKF on a nonlinear benchmark for system identification, the cascaded tanks data set. 
A detailed description can be found in <cit.>.
The training and testing data set consists of 1024 data points.
To train our GP model, we choose lagged inputs and outputs as input to our GP, as described in <cit.>, resulting in an input of dimensionality D=14.
We use a squared exponential kernel, which hyperparameters we optimize with the Gaussian process toolbox by <cit.>, and we choose I=4, such that the model has M=4^14=268435456 parameters.

For the comparison to the TNKF, we choose the TT-ranks for the mean to be R_2=R_14=4, R_3=⋯ R_13=10, and we vary R_ and the TTm-ranks of the covariance matrix for the TNKF denoted by R_.
Fig. <ref> and <ref> show the RMSE and NLL over the time iterations of the respective filter. 
When R_=1 and R_=1, both methods perform almost the same, as visualized by the overlapping green and blue lines.
When R_^2=R_=4, our method improves both prediction accuracy and uncertainty estimation compared to the R_=1.
On the contrary, the TNKF diverges and leaves the plotted figure area because the covariance matrix loses positive definiteness.
When R_^2=R_=16, the TNKF shows a similar behavior, while the TNSRKF results in lower RMSEs but mostly higher NLL values.
This setting shows that higher values for R_ are not always beneficial for the uncertainty estimation.

Finally, Fig. <ref> shows the predictions with the TNSRKF on testing data after seeing 100, 200, and 922 data points.
Aligned with the plot showing the RMSE and NLL, after 100 data points, the prediction is quite bad and uncertain.
After 200 data points, the predictions are better and more certain and further improve after seeing the entire data set.









§ CONCLUSION

In this paper, we presented a TT-based solution for online GP regression in terms of an SRKF.
In our experiments, we show that our method is scalable to a high number of input dimensions at a reasonable computational cost such that all experiments could be run on a conventional laptop.
In addition, we improve the state-of-the-art method for TN-based Kalman filter: In settings where the TNKF loses positive (semi-)definiteness and becomes numerically unstable, our method avoids this issue because we compute the square root covariance factors instead of the covariance matrix.
In this way, we can choose settings for our method that achieve better accuracy than the TNKF.

A future work direction is online hyperparameter optimization.
We are looking at a truly online scenario, so future data is not available.
Thus we cannot swipe over mini-batches of data multiple times like other methods, e.g. <cit.>, to optimize hyperparameters.

Finally, there is still ongoing research to determine how to choose TT-ranks and TTm-ranks.
In the synthetic experiments, we showed the impact of R_ and R_.
Generally, the TT- and TTm-ranks need to be treated as hyperparameters.




plain

