METHOD FOR CALIBRATING ON-LINE AND WITH FORGETTING FACTOR A DIRECT NEURAL INTERFACE WITH PENALISED MULTIVARIATE REGRESSION

Abstract

The present invention relates to a method for calibrating on-line a direct neural interface implementing a REW-NPLS regression between an output calibration tensor and an input calibration tensor. The REW-NPLS regression comprises a PARAFAC iterative decomposition of the cross covariance tensor between the input calibration tensor and the output calibration tensor, each PARAFAC iteration comprising a sequence of M elementary steps (2401, 2401, . . . 240M) of minimisation of a metric according to the alternating least squares method, each elementary minimisation step relating to a projector and considering the others as constant, said metric comprising a penalisation term that is a function of the norm of this projector, the elements of this projector not being subjected to a penalisation during a PARAFAC iteration f not being penalisable during following PARAFAC iterations. Said calibration method makes it possible to obtain a predictive model of which the non-zero coefficients are sparse blockwise.

Description

TECHNICAL FIELD

The present invention relates to the field of direct neural interfaces, also called BCI (Brain Computer Interfaces) or BMI (Brain Machine Interfaces). It is notably applied to the direct neural command of a machine, such as an exoskeleton or a computer.

PRIOR ART

Direct neural interfaces use the electrophysiological signals emitted by the cerebral cortex to elaborate a command signal. These neural interfaces have been the subject of numerous research notably with the aim of restoring a motor function to a paraplegic or tetraplegic subject using a motorised prothesis or orthosis.

Neural interfaces may be of invasive or non-invasive nature. Invasive neural interfaces use intracortical electrodes (that is to say implanted in the cortex) or cortical electrodes (arranged on the surface of the cortex) collecting in this latter case electrocorticographic (EcoG) signals. Non-invasive neural interfaces use electrodes placed on the scalp to collect electrocephalographic (EEG) signals. Other types of sensors have also been envisaged such as magnetic sensors measuring the magnetic fields induced by the electrical activity of the neurones of the brain. One then speaks of magnetocephalographic (MEG) signals.

Direct neural interfaces advantageously use ECoG type signals, having the advantage of a good compromise between biocompatibility (matrix of electrodes implanted on the surface of the cortex) and quality of the collected signals.

The ECoG signals thus measured must be processed in order to estimate the trajectory of the movement desired by the subject and to deduce therefrom the command signals of the computer or of the machine. For example, when it involves commanding an exoskeleton, the BCI estimates the trajectory of the desired movement from measured electrophysiological signals and deduces therefrom the control signals enabling the exoskeleton to reproduce the trajectory in question. Similarly, when it involves commanding a computer, the BCI estimates for example the desired trajectory of a pointer or a cursor from electrophysiological signals and deduces therefrom the command signals of the cursor/pointer.

The estimation of trajectory, and more specifically that of kinematic parameters (position, speed, acceleration), is also designated neural decoding in the literature. Neural decoding notably makes it possible to command a movement (of a prothesis or a cursor) from ECoG signals.

The estimation of trajectory and the calculation of signals for controlling the exoskeleton or the effector generally requires a prior learning or calibration phase, designated off-line. During this phase, the subject imagines, observes or makes a movement according to a determined trajectory during a given calibration interval. The electrophysiological signals measured during this interval are exploited in relation with this trajectory to construct a predictive model and more specifically to calculate the parameters of this model.

The validity of the predictive model is however limited in time on account of the non-stationarity of the neural signals. For this reason, it is necessary to carry out an on-line calibration of the predictive model, that is to say as and when the neural signals are observed and the command applied.

Different methods for calibrating on-line a neural interface have been described in the prior art. However, these on-line calibration methods may require the manipulation of a considerable amount of data resulting from the juxtaposition of data from prior calibration steps and those from a new calibration step.

An on-line BCI calibration method has been described in the article of A. Eliseyev et al. entitled “Recursive exponentially weighted N-way Partial Least Squares regression with recursive validation of hyper-parameters in Brain-Computer Interface applications” published in Scientific Reports, vol. 7, no 1, p. 16281, November 2017 as well as in the patent application FR-A-3 061 318. This method will be designated hereafter under the acronym REW-NPLS (Recursive Exponentially Weighted N-way Partial Least Squares).

Furthermore, different alternatives of off-line calibration method have been proposed in the literature, with the aim of improving the robustness of the prediction and to increase the sparsity of non-zero coefficients in the predictive model tensor. In particular, the sparsity of non-zero coefficients in the spatial mode of this tensor makes it possible to only select certain relevant electrodes for the prediction of the trajectory and the calculation of the command. Such an alternative using a penalised multipath regression with such an aim has been described in the article of A. Eliseyev et al. entitled “L1-penalised N-way PLS for subset of electrodes selection in BCI experiments” published in J. Neural Eng. vol. 9, no 4, p. 045010 as well as in the patent application FR-A-3 046 471.

However, such a calibration method using a penalised multipath regression cannot be implemented on-line on account of the considerable amount of calculations that it requires.

The subject matter of the present invention is consequently to propose a method for calibrating on-line a direct neural interface that can also increase the sparsity of non-zero coefficients of the prediction tensor, notably in such a way as to only select at each observation window a relevant sub-set of electrodes.

DESCRIPTION OF THE INVENTION

The present invention is defined by a method for calibrating on-line a direct neural interface intended to receive a plurality of electrophysiological signals during a plurality of observation windows associated with observation times, to form an input tensor and to deduce therefrom, by means of a predictive model, an output tensor providing command signals intended to command one or more effectors, said calibration is carried out during a plurality of calibration steps, each calibration step u being carried out from input calibration data, represented by an input calibration tensor, X_u∈□^{(I1× . . . ×IN)}, and from output calibration data represented by an output calibration tensor Y_u∈□^{(J1× . . . ×JM)}, said calibration step implementing a REW-NPLS regression comprising:

the calculation of a covariance tensor XX_u from a covariance tensor XX_u−1 obtained during a preceding calibration step and the tensor X_u, as well as the calculation of a cross covariance tensor XY_u from a cross covariance tensor XY_u−1 obtained during said preceding calibration step and the tensors X_u and Y_u;

a PARAFAC iterative decomposition of the cross covariance tensor XY_u, each PARAFAC iteration being associated with a dimension f of the latent variables space on which said regression is carried out, each iteration providing a plurality of projectors w_fⁱ, i=1, . . . , M of respective sizes I₁, I₂, . . . I_M and a predictive model defined by a prediction coefficients tensor B_u^f∈□^{(I1× . . . ×IM)×(J1× . . . ×JM)} and a bias tensor β_u^f∈□^{(J1× . . . ×JM)}.

Said calibration method is original in that each PARAFAC iteration comprises a sequence of M elementary steps of minimisation of a metric according to the alternating least squares method, each elementary minimisation step relating to a projector and considering the others as constant, said metric comprising a penalisation term that is a function of the norm of this projector, the elements of this projector not being subjected to a penalisation during a PARAFAC iteration f not being penalisable during following PARAFAC iterations.

Advantageously, an elementary minimisation step relating to a projector w_fⁱ aims to minimise the penalised metric

$\underset{{\tilde{w}}_f^i}{\min }\left({{\underset{\_}{\mathrm{XY}}}_{u\left(i\right)}-{{\tilde{w}}_f^i\left(w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1\right)}^T}^2+{\lambda }_i{{\tilde{w}}_f^i}_{q,{\Omega }_{i,f}}\right)$

where {tilde over (w)}_fⁱ is a non-normalised projector and

$w_f^m=\frac{{\tilde{w}}_f^m}{{\tilde{w}}_f^m},$

m=1, . . . ; m≠i are normalised projectors, XY_u(i) is a matrix of size

$I_i\times \prod _{\underset{n\ne i}{n=1}}^N{I}_n$

obtained by deployment of the cross covariance tensor XY_u while taking the mode i as reference, ⊗ is the Kronecker product, λ_i is a strictly positive coefficient, ∥w∥_q,Ωi,f, q=0, ½, 1 is a norm defined by

${w_f^i}_{0,{\Omega }_{i,f},}=\sum _{k\in {\Omega }_{i,f}}\left(1-{\delta }_{0,w_{f,k}^i}\right),{w_f^i}_{1/2,{\Omega }_{i,f}}=\sum _{k\in {\Omega }_{i,f}}\sqrt{w_{f,k}^i},{w_f^i}_{1,{\Omega }_{i,f}}=\sum _{k\in {\Omega }_{i,f}}w_{f,k}^i$

where w_f,kⁱ is the k^th element of the projector w_fⁱ, δ is the Kronecker symbol, and Ω_i,f is the set of indices corresponding to the penalisable elements of the projector w_fⁱ at the PARAFAC iteration f.

Alternatively, an elementary minimisation step relating to a projector w_fⁱ aims to minimise the penalised metrics

$\underset{w_{f,k}^i}{\min }\left({z_i^k-{w_{f,k}^i\left({\left(w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1\right)}^T\right)}^T}^2+{\lambda }_i{g}_q\left(w_{f,k}^i\right)\right),$

k=1, . . . , I_i where z_i^k is the k^th column vector of the matrix XY_u(i) of size

$I_i\times \prod _{\underset{n\ne i}{n=1}}^N{I}_n$

obtained by deployment of the cross covariance tensor XY_u while taking the mode i as reference, ⊗ is the Kronecker product, λ_i is a strictly positive coefficient, and the functions g_q, q=0, ½, 1 are defined by:

g ₀(w_f,kⁱ)=1−δ_0,wf,ki if k∈Ω_i,f and g₀(w_f,kⁱ)=0 if not

g _1/2(w_f,kⁱ)=√{square root over (|w_f,kⁱ|)} if k∈Ω_i,f and g_1/2(w_f,kⁱ)=0 if not

g ₁(w_f,kⁱ)=|w_f,kⁱ| if k∈Ω_i,f and g₁(w_f,kⁱ)=0 if not

where Ω_i,f is the set of indices k corresponding to the penalisable elements w_f,kⁱ of the projector w_fⁱ at the PARAFAC iteration f.

In a particular embodiment, q=0 and the elements w_f,kⁱ resulting from the minimisation are given by:

(w_f,kⁱ)_L0=0 if k∈Ω_i,f and |(w_f,kⁱ)_LS|≤Th_L0ⁱ

(w_f,kⁱ)_L0=(w_f,kⁱ)_LS if not

where (w_f,kⁱ)_LS represents the minimisation solution in the least squares sense:

${\left(w_{f,k}^i\right)}_{LS}=\frac{z_i^k\left(w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1\right)}{{w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1}^2}$

and that the threshold Th_L0ⁱ is given by

${\mathrm{Th}}_{L0}^i=\frac{\sqrt{{\lambda }_i}}{w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1}.$

In a second particular embodiment, q=½ and the elements w_f,kⁱ resulting from the minimisation are given by:

(w_f,kⁱ)_L0=0 if k∈Ω_i,f and |(w_f,kⁱ)_LS|≤Th_L0,5ⁱ

(w_f,kⁱ)_L0=arg min(h_0,5(0),h_0,5(Bⁱ·(w_f,kⁱ)_LS)) if k∈Ω_i,f and |(w_f,kⁱ)_LS|>Th_L0,5ⁱ

(w_f,kⁱ)_L0=(w_f,kⁱ)_LS if not

where h_0,5(w_f,kⁱ)=∥z₁^k−w_f,kⁱ(w_f³≃w_f²)^T∥²+λ_i√{square root over (|w_f,kⁱ|)} and where (w_f,kⁱ)_LS represents the minimisation solution in the least squares sense:

${\left(w_{f,k}^i\right)}_{LS}=\frac{z_i^k\left(w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1\right)}{{w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1}^2}$

and the threshold Th_L0,5ⁱ is given by

${\mathrm{Th}}_{L0,5}^i=\frac{3}{4}{\left(\frac{{\lambda }_i}{{w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1}^2}\right)}^{\frac{2}{3}}.$

In a third particular embodiment, q=1 and the elements w_f,kⁱ resulting from the minimisation are given by:

(w_f,kⁱ)_L1=0 if k∈Ω_i,f and |(w_f,kⁱ)_LS|≤Th_L1ⁱ

(w_f,kⁱ)_L1=sgn((w_f,kⁱ)_LS)(|(w_f,kⁱ)_LS|−Th_L1ⁱ) if k∈Ω_i,f and |(w_f,kⁱ)_LS|>Th_L1ⁱ

(w_f,kⁱ)_L0=(w_f,kⁱ)_LS if not

where (w_f,kⁱ)_LS represents the minimisation solution in the least squares sense:

${\left(w_{f,k}^i\right)}_{LS}=\frac{z_i^k\left(w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1\right)}{{w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1}^2}$

and the threshold Th_L1ⁱ is given by

${\mathrm{Th}}_{L1}^i=\frac{{\lambda }_i}{w_f^M\otimes \dots w_f^{i+1}\otimes w_f^{i-1}\otimes \dots \otimes w_f^1}$

Typically, said on-line calibration method may apply to electrophysiological signals constituted of ECoG signals.

BRIEF DESCRIPTION OF THE DRAWINGS

Other characteristics and advantages of the invention will become clear on reading the preferential embodiment of the invention, described with reference to the appended figures among which:

FIG. 1A schematically represents the decomposition of a tensor by parallel factorial analysis (PARAFAC);

FIG. 1B schematically represents the development of a tensor with respect to a mode;

FIG. 1C schematically represents a decomposition using an alternating least squares (ALS) method;

FIG. 2 schematically represents the flowchart of a recursive method for calibrating a direct neural interface with penalised multivariate regression, according to an embodiment of the present invention;

FIG. 3 schematically represents the mechanism for updating the predictive model in the calibration method of FIG. 2.

DETAILED DESCRIPTION OF PARTICULAR EMBODIMENTS

A direct neural interface with continuous decoding will be considered hereafter.

The electrophysiological signals coming from the different electrodes are sampled and grouped together by data blocks, each block corresponding to a sliding observation window of width ΔT. Each observation window is defined by an observation time (epoch) at which the window in question starts.

The electrophysiological signals may be subjected to a pre-processing. This pre-processing may notably comprise a deletion of the average taken over the set of electrodes, then a time-frequency analysis is carried out on each of the observation windows.

The time-frequency analysis is based on a decomposition into wavelets, notably Morlet wavelets. Those skilled in the art will understand however that other types of time-frequency analysis could be envisaged by those skilled in the art.

These time-frequency analyses results may further be subjected to frequency smoothing or decimation.

Thus, with each observation window, or observation time t, is associated a tensor of order 3 of observation data, hence the generation of an input tensor of order 4: the first mode corresponds to the successive observation windows, the second mode corresponds to the space, in other words to the sensors, the third mode corresponds to the time within an observation window, in other words to the positions of the wavelets, and the fourth mode corresponds to the frequency, in other words to the number of frequency bands used for the decomposition into wavelets over an observation window.

More generally, the input tensor (or observation tensor) will be of order N, the first mode being in all cases that relative to the observation times (epochs). The input tensor (or observation tensor) is noted X and is of dimension I₁× . . . ×I_N.

In the same way, the trajectory of the movement imagined, observed or made is described by an output tensor (or command tensor) of order M, noted Y, of dimension J₁× . . . ×J_M, of which the first mode corresponds to the successive times to which the commands will apply (as a general rule, this first mode also corresponds to the observation windows), the other modes corresponding to the commands of different effectors or to the other degrees of freedom of a multi-axis robot.

More specifically, the output tensor provides command data blocks, each of the blocks making it possible to generate command signals relative to the different effectors or degrees of freedom. Thus, it will be understood that the dimension of each data block could depend on the envisaged case of use and notably the number of degrees of freedom of the effector.

One will note hereafter X, the observation tensor at the time t. This tensor is consequently of order N and of dimension I₁× . . . ×I_N. It takes its values in a space X⊂□^{I1× . . . ×IN} where □ is the set of real values. Similarly, one will note Y, the command tensor at the time t. This output tensor is of order M and of dimension J₁× . . . ×J_M. It takes its values in a space Y⊂□X_u∈□^{J1× . . . ×JM}.

The calibration or, in an equivalent manner, the mechanism for updating the predictive model according to the REW-NPLS (Recursive Exponentially Weighted N-way Partial Least Squares) algorithm is recalled hereafter.

If one notes {X_u, Y_u} the calibration data blocks (or training data sets) during the iteration u the basic idea of the REW-NPLS algorithm is to update the predictive model from the covariance and cross covariance tensors obtained at the preceding calibration step, respectively noted XX_u−1 and XY_u−1. These two tensors make it possible to obtain a condensed representation of the predictive model at the preceding step, without there being a need to store the history of preceding calibration data blocks: {X₁, Y₁}, {X₂, Y₂}, . . . {X_u−1, Y_u−1}.

More specifically, the covariance and cross covariance tensors are updated by means of recurrence relationships:

XX _u =γXX _u−1 +X _u×₁X_u  (1-1)

XY _u =γXY _u−1 +X _u×₁Y_u  (1-2)

where ×₁ designates the tensor product of mode 1 and γ is a forgetting coefficient with 0<γ<1.

One thus understands that a current calibration period, carried out during the iteration u takes account of a preceding calibration period carried out during the iteration u−1 and the tensors of the observation data of the current calibration period, the history being weighted by the forgetting coefficient γ with 0<γ<1.

This operation of updating the covariance tensors is carried out in an iterative manner over the calibration periods.

For each update, the predictive model is next obtained in an iterative manner by means of a PARAFAC decomposition (generalisation to the tensorial case of a decomposition into singular values) of the cross covariance tensor, XY_u. It will be recalled that, generally speaking, a PARAFAC decomposition of a tensor Z of order M, i.e. Z∈□^{I1×I2× . . . ×IM} makes it possible to represent this tensor in the form of a linear combination of external products of vectors (tensors of order 1) to a nearest residual tensor:

$\begin{array}{cc}\underset{\_}{Z}=\sum _{r=1}^R{\theta }_r{w}_r^1\circ w_r^2\circ \dots \circ w_r^M+\underset{\_}{E}& \left(2\right)\end{array}$

with ∥w_r^m∥=1, m=1, . . . , M where ∘ is the external product and E is a residual tensor and _r is the index of the PARAFAC iteration and R is a given integer. A PARFAC decomposition of a tensor Z has been represented schematically in FIG. 1A.

In the case of a REW-NPLS calibration step, each PARAFAC iteration only relates to an elementary decomposition, that is to say R=1. It makes it possible to extract from the tensor XY_u a set of projection vectors, also designated projectors, w_f¹, . . . , w_f^M, of respective I₁, I₂, . . . , I_M where the iteration number f gives the dimension of a latent variables space onto which the input tensor X_u is projected.

In order to simplify the presentation and without prejudice of generalisation to a tensor of any order, we will assume hereafter that the order of the tensor XY_u is equal to 3 (M=3), that is to say XY_u∈□^I1×I2×I3 with ∥XY_u∥=1 For example, the first mode corresponds to the epoch, the second to the frequency and the third to the space (that is to say to the electrodes).

At the iteration f of the PARAFAC decomposition, one searches for the projection vectors w_f¹, w_f², w_f³ which verify:

$\underset{}{\min }{{\underset{\_}{\mathrm{XY}}}_u-{\underset{\_}{Y}}_u}^2\mathrm{with}$ ${\underset{\_}{Y}}_u={\theta }_f{w}_f^1\circ w_f^2\circ w_f^3\mathrm{and}$ $w_f^1=w_f^2=w_f^3=1$

This minimisation problem may be resolved by an alternating least squares method in which one carries out sequentially a minimisation for each of the vectors of the decomposition, i.e.:

$\begin{array}{cc}\underset{w_f^1}{\min }{{\underset{\_}{\mathrm{XY}}}_{u\left(1\right)}-{w_f^1\left(w_f^3\otimes w_f^2\right)}^T}^2& \left(4\text{-}1\right)\\ \underset{w_f^2}{\min }{{\underset{\_}{\mathrm{XY}}}_{u\left(2\right)}-{w_f^2\left(w_f^3\otimes w_f^1\right)}^T}^2& \left(4\text{-}2\right)\\ \underset{w_f^3}{\min }{{\underset{\_}{\mathrm{XY}}}_{u\left(3\right)}-{w_f^3\left(w_f^2\otimes w_f^1\right)}^T}^2& \left(4\text{-}3\right)\end{array}$

where ⊗ is the Kronecker product and Z_(n) is a matrix obtained by deployment of the tensor Z while taking the mode n as reference mode (tensor unfolding or flattening along mode n). For example, the matrix Z₍₁₎ is defined by (z₁¹| . . . |z₁^I1)∈□^I1×I2I3 where the Z₁^j are the columns of the matrix Z₍₁₎.

The minimisation operations (4-1), (4-2) and (4-3) are repeated up to convergence. the minimisation solutions in the least squares (LS) sense for each of these operations are given by the projectors:

$\begin{array}{cc}{w}_f^1=\frac{{\underset{\_}{\mathrm{XY}}}_{u\left(1\right)}\left(w_f^3\otimes w_f^2\right)}{{w_f^3\otimes w_f^2}^2}& \left(5\text{-}1\right)\\ w_f^2=\frac{{\underset{\_}{\mathrm{XY}}}_{u\left(2\right)}\left(w_f^3\otimes w_f^1\right)}{{w_f^3\otimes w_f^1}^2}& \left(5\text{-}2\right)\\ w_f^3=\frac{{\underset{\_}{\mathrm{XY}}}_{u\left(3\right)}\left(w_f^2\otimes w_f^1\right)}{{w_f^2\otimes w_f^1}^2}& \left(5\text{-}3\right)\end{array}$

Each iteration f of the PARAFAC decomposition provides a plurality M of projectors w_f¹, . . . , w_f^M and a prediction coefficients tensor, or predictive model of the direct neural interface, defined by a prediction coefficients tensor B_N^f∈□^{(I1× . . . ×IN)×(J1× . . . ×JM)} and a bias tensor β_u^f∈^{(J1× . . . ×JM)}, making it possible to predict the command tensor from the observation tensor. One thus obtains a plurality of models, B_u^f, f=1, . . . , F corresponding to different successive dimensions of the latent variables space.

One next determines, at the following calibration step, u+1, among B_u^f, f=1, . . . , F, the predictive model B_u^f* corresponding to the smallest prediction error where one has noted f* the index f corresponding to this smallest prediction error. This determination of the predictive model at the following calibration step is designated in the literature under the term recursive validation.

FIG. 1C schematically represents a decomposition using an alternating least squares (ALS) method.

The decomposition by alternating least squares has been illustrated for a REW-NPLS calibration step and a PARAFAC iteration f.

At step 110, one estimates the vector w¹_f (corresponding to the epoch mode) from the matrix XY_u(1) and the current vectors w_f² and w_f³. This operation may be considered as a projection of the column vectors of the matrix XY_u(1) onto the tensorial space generated by the modes 2 and 3.

At step 120, one estimates the vector w²_f (corresponding to the frequency mode) from the matrix XY_u(2) and the current vectors w_f¹ and w_f³. This operation may be considered as a projection of the column vectors of the matrix XY_u(2) onto the tensorial space generated by the modes 1 and 3.

At step 130, one estimates the vector w³_f (corresponding to the spatial mode) from the matrix XY_u(3) and the current vectors w_f¹ and w_f². This operation may be considered as a projection of the column vectors of the matrix XY_u(3) onto the tensorial space generated by the modes 1 and 2.

Steps 110, 120, 130 are iterated until a convergence criterion is met, for example when the sum of the deviations between the vectors w_f^j of two successive ALS iterations is less than a predetermined threshold value.

The basic idea of the present invention is to carry out a penalisation by section at each iteration of the PARAFAC algorithm, a penalisation by section resulting in an elimination of elements relative to an index value in the tensor of the prediction model. The elements not being subjected to a penalisation during an iteration f are no longer penalisable for the following iterations, f+1, f+2, . . . , F. On the other hand, the elements penalised during an iteration are taken into account during the PARAFAC decomposition of the following iteration. One thus promotes a blockwise sparsity of non-zero elements in the tensor of the prediction model, which considerably simplifies the calculation of the command.

More specifically, during an iteration f of the PARAFAC decomposition, the minimisation of the quadratic deviations according to (4-1), (4-2) and (4-3) in the alternating least squares method is replaced by a minimisation taking into account a penalisation term favouring this sparsity, namely:

$\begin{array}{cc}\underset{{\tilde{w}}_f^1}{\min }({\underset{\_}{\mathrm{XY}}}_{u\left(1\right)}-{{\tilde{w}}_f^1\left(w_f^3\otimes w_f^2\right)}^T{}^2+{\lambda }_1{{\tilde{w}}_f^1}_{q,{\Omega }_{1,f}})& \left(6\text{-}1\right)\\ \underset{{\tilde{w}}_f^2}{\min }({\underset{\_}{\mathrm{XY}}}_{u\left(2\right)}-{{\tilde{w}}_f^2\left(w_f^3\otimes w_f^1\right)}^T{}^2+{\lambda }_2{{\tilde{w}}_f^2}_{q,{\Omega }_{2,f}})& \left(6\text{-}2\right)\\ \underset{{\tilde{w}}_f^3}{\min }({\underset{\_}{\mathrm{XY}}}_{u\left(3\right)}-{{\tilde{w}}_f^3\left(w_f^2\otimes w_f^1\right)}^T{}^2+{\lambda }_3{{\tilde{w}}_f^3}_{q,{\Omega }_{3,f}})& \left(6\text{-}3\right)\end{array}$

where {tilde over (w)}_f¹, {tilde over (w)}_f² and {tilde over (w)}_f³ are non-normalised projectors and

$w_f^1=\frac{{\tilde{w}}_f^1}{{\tilde{w}}_f^1},w_f^2=\frac{{\tilde{w}}_f^2}{{\tilde{w}}_f^2},w_f^3=\frac{{\tilde{w}}_f^3}{{\tilde{w}}_f^3}$

are the corresponding normalised projectors and λ₁, λ₂, λ₃ are coefficients that are strictly positive and less than 1. It will be noted that due to the presence of penalisation terms in the expressions (6-1), (6-2) and (6-3), the minimisation can no longer be carried out on normalised projection vectors as in the classic ALS method. The norms ∥w_fⁱ∥_q,Ωi q=0, ½, 1 are defined by:

$\begin{array}{cc}{w_f^i}_{0,{\Omega }_{i,f},}=\sum _{k\in {\Omega }_{i,f}}\left(1-{\delta }_{0,w_{f,k}^i}\right)& \left(7\text{-}1\right)\\ {w_f^i}_{1/2,{\Omega }_{i,f}}=\sum _{k\in {\Omega }_{i,f}}\sqrt{w_{f,k}^i}& \left(7\text{-}2\right)\\ {w_f^i}_{1,{\Omega }_{i,f}}=\sum _{k\in {\Omega }_{i,f}}w_{f,k}^i& \left(7\text{-}3\right)\end{array}$

where w_f,kⁱ is the k^th element of the projection vector w_fⁱ, δ is the Kronecker symbol, and Ω_i,f is the set of indices corresponding to the penalisable elements, it being understood that Ω_i,f+1∩Ω_i,f since the non-penalised elements of a projector at an iteration f are not penalisable at the iteration f+1.

It will notably be understood that the norm of the expression (7-1) gives the sum of the non-zero elements of which the indices belong to Thus for example, in this case, the penalisation will be all the greater when the number of elements of which the indices belong to Ω_i,j is higher.

In order to simplify the minimisation calculation in the expressions (6-1), (6-2), (6-3), one could advantageously operate element by element, namely:

$\begin{array}{cc}\underset{w_{f,k}^1}{\min }\left({z_1^k-{w_{f,k}^1\left(w_f^3\otimes w_f^2\right)}^T}^2+{\lambda }_1{g}_q\left(w_{f,k}^1\right)\right),k=1,\dots ,I_1& \left(8\text{-}1\right)\\ \underset{w_{f,k}^2}{\min }\left({z_2^k-{w_{f,k}^2\left(w_f^3\otimes w_f^1\right)}^T}^2+{\lambda }_2{g}_q\left(w_{f,k}^2\right)\right),k=1,\dots ,I_2& \left(8\text{-}2\right)\\ \underset{w_{f,k}^3}{\min }\left({z_3^k-{w_{f,k}^3\left(w_f^2\otimes w_f^1\right)}^T}^2+{\lambda }_3{g}_q\left(w_{f,k}^3\right)\right),k=1,\dots ,I_3& \left(8\text{-}3\right)\end{array}$

where the vectors z₁^k, z₂^k, z₃^k are the k^th column vectors of the matrices XY_u(1), XY_u(2), XY_u(3) obtained by deployment of the cross covariance tensor XY_u according to modes 1, 2 and 3, respectively.

The functions g_q, q=0, ½, 1 apply to scalar values and are defined by:

g ₀(w_f,kⁱ)=1−δ_0,wf,ki if k∈Ω_i,f and g₀(w_f,kⁱ)=0 if not  (9-1)

g _1/2(w_f,kⁱ)=√{square root over (|w_f,kⁱ|)} if k∈Ω_i,f and g_1/2(w_f,kⁱ)=0 if not  (9-2)

g ₁(w_f,kⁱ)=|w_f,kⁱ| if k∈Ω_i,f and g₁(w_f,kⁱ)=0 if not  (9-3)

It may be shown that the result of the penalised minimisation operation according to (8-1), (8-2), (8-3) for the norm L₀ (that is to say q=0) gives:

(w_f,kⁱ)_L0=0 if k∈Ω_i,f and |(w_f,kⁱ)_LS|≤Th_L0ⁱ  (10-1)

(w_f,kⁱ)_L0=(w_f,kⁱ)_LS if not  (10-2)

where (w_f,kⁱ)_LS represents the minimisation solution in the least squares sense (5-1), (5-2), (5-3), that is to say:

$\begin{array}{cc}{\left(w_{f,k}^1\right)}_{LS}=\frac{z_1^k\left(w_f^3\otimes w_f^2\right)}{{w_f^3\otimes w_f^2}^2}& \left(11\text{-}1\right)\\ {\left(w_{f,k}^2\right)}_{LS}=\frac{z_2^k\left(w_f^3\otimes w_f^1\right)}{{w_f^3\otimes w_f^1}^2}& \left(11\text{-}2\right)\\ {\left(w_{f,k}^3\right)}_{LS}=\frac{z_3^k\left(w_f^2\otimes w_f^1\right)}{{w_f^2\otimes w_f^1}^2}& \left(11\text{-}3\right)\end{array}$

and the thresholds Th_L0ⁱ, i=1, 2, 3 are given by:

$\begin{array}{cc}T{h}_{L0}^1=\frac{\sqrt{{\lambda }_1}}{w_f^3\otimes w_f^2}& \left(12\text{-}1\right)\\ T{h}_{L0}^2=\frac{\sqrt{{\lambda }_2}}{w_f^3\otimes w_f^1}& \left(12\text{-}2\right)\\ T{h}_{L0}^3=\frac{\sqrt{{\lambda }_3}}{w_f^2\otimes w_f^1}& \left(12\text{-}3\right)\end{array}$

The penalised minimisation operation consequently comes down to carrying out a thresholding on the elements of the projectors obtained at each step of the ALS method.

Similarly, it may be shown that the penalised minimisation operation according to (8-1), (8-2), (8-3) for the norm L_1/2 (that is to say q=½) has for solution:

(w_f,kⁱ)_L0=0 if k∈Ω_i,f and |(w_f,kⁱ)_LS|≤Th_L0,5ⁱ  (13-1)

(w_f,kⁱ)_L0=arg min(h_0,5(0),h_0,5(Bⁱ·(w_f,kⁱ)_LS)) if k∈Ω_i,f and |(w_f,kⁱ)_LS|>Th_L0,5ⁱ   (13-2)

(w_f,kⁱ)_L0=(w_f,kⁱ)_LS if not  (13-3)

where h_0,5(w_f,kⁱ)=∥z₁^k−w_f,kⁱ(w_f³≃w_f²)^T∥²+λ_i√{square root over (|w_f,kⁱ|)} and the thresholds Th_L0,5ⁱ, i=1, 2, 3 are given by:

$\begin{array}{cc}{\mathrm{Th}}_{L0,5}^1=\frac{3}{4}{\left(\frac{{\lambda }_1}{{w_f^3\otimes w_f^2}^2}\right)}^{\frac{2}{3}}& \left(14\text{-}1\right)\\ {\mathrm{Th}}_{L0,5}^2=\frac{3}{4}{\left(\frac{{\lambda }_2}{{w_f^3\otimes w_f^1}^2}\right)}^{\frac{2}{3}}& \left(14\text{-}2\right)\\ {\mathrm{Th}}_{L0,5}^3=\frac{3}{4}{\left(\frac{{\lambda }_3}{{w_f^2\otimes w_f^1}^2}\right)}^{\frac{2}{3}}& \left(14\text{-}3\right)\end{array}$

The coefficients Bⁱ, i=1, 2, 3 are defined as the respective solutions of the cubics x(1−x)²=C_i, i=1, 2, 3, with:

$\begin{array}{cc}{C}^1=\frac{{\lambda }_1^2}{16{w_f^3\otimes w_f^2}^4{\left({\left(w_{f,k}^1\right)}_{\mathrm{LS}}\right)}^3}& \left(15\text{-}1\right)\\ C^2=\frac{{\lambda }_2^2}{16{w_f^3\otimes w_f^1}^4\left({\left(w_{f,k}^2\right)}_{\mathrm{LS}}\right)}& \left(15\text{-}2\right)\\ C^3=\frac{{\lambda }_3^2}{16{w_f^2\otimes w_f^1}^4{\left({\left(w_{f,k}^3\right)}_{\mathrm{LS}}\right)}^3}& \left(15\text{-}3\right)\end{array}$

Finally, similarly, it may be shown that the penalised minimisation operation according to (8-1), (8-2), (8-3) for the norm L₁ (that is to say q=1) has for solution:

(w_f,kⁱ)_L1=0 if k∈Ω_i,f and |(w_f,kⁱ)_LS|≤Th_L1ⁱ   (16-1)

(w_f,kⁱ)_L1=sgn((w_f,kⁱ)_LS)(|(w_f,kⁱ)_LS|−Th_L1ⁱ) if k∈Ω_i,f and |(w_f,kⁱ)_LS|>Th_L1ⁱ   (16-2)

(w_f,kⁱ)_L0=(w_f,kⁱ)_LS if not  (16-3)

where the thresholds Th_L1ⁱ, i=1, 2, 3 are given by:

$\begin{array}{cc}{\mathrm{Th}}_{L1}^1=\frac{{\lambda }_1}{w_f^3\otimes w_f^2}& \left(17\text{-}1\right)\\ {\mathrm{Th}}_{L1}^2=\frac{{\lambda }_2}{w_f^3\otimes w_f^1}& \left(17\text{-}2\right)\\ {\mathrm{Th}}_{L1}^3=\frac{{\lambda }_3}{w_f^2\otimes w_f^1}& \left(17\text{-}3\right)\end{array}$

It will be understood that, whatever the norm used, the penalised minimisation is advantageously reduced to a simple thresholding of the elements of a projector at each step of the ALS method. This penalised minimisation makes it possible to promote blockwise sparsity of non-zero coefficients in the predictive model and consequently to simplify the calculation. It is important to note that this may be carried out on-line with the aim of carrying out an incremental calibration of the predictive model.

FIG. 2 schematically represents the flow chart of a recursive method for calibrating a direct neural interface with penalised multivariate regression, according to an embodiment of the present invention.

At step 200, one takes into account new input and output calibration data, respectively in the form of an input tensor X_u and an output tensor Y_u. The input calibration data result from the measurement of electrophysiological signals. The output calibration data result for example from the measurement of kinematic parameters of the trajectory or command signals of effectors. The input and output tensors may be if appropriate subjected to a pre-processing aiming to centre and normalise them, in a manner known per se.

At step 210, one determines among the set of predictive models B_u−1^f, β_u−1^f, f=1, . . . , F obtained at the preceding step, the predictive model {B_u−1^f*, β_u−1^f*} corresponding to the lowest prediction error. This model is used to calculate the command tensor (output tensor) from the observation tensor (input tensor).

At step 220, one updates the covariance XX_u and cross covariance, XY_u, tensors according to the relationships (1-1) and (1-2), while taking account of the forgetting factor.

One next enters into a first NPLS iterative loop, in which one iterates on the dimension f of the latent variables space, with f=1, . . . , F where F is a maximum dimension. This iterative loop makes it possible to generate successively a plurality of predictive models of the direct neural interface, each predictive model being defined by a prediction coefficients tensor B_u^f∈□^{(I1× . . . ×IM)×(J1× . . . ×JM)} and a bias tensor β_u^f∈□^{(f1× . . . ×JM)}.

At step 230, one initialises f with f=1 and the sets Ω_i,1, i=1, . . . , M with Ω_i,1={1, 2, . . . , I_i} for i=1, . . . , M.

For each value of f, one carries out a PARAFAC decomposition of the cross covariance tensor using an alternating least squares method. Unlike a conventional PARAFAC decomposition (cf. expressions (4-1), (4-2), (4-3)) used in the REW-NPLS method, a PARAFAC decomposition regularised by means of a penalised metric is used here, the penalisation term aiming to promote blockwise sparsity of non-zero coefficients within the tensor representative of the predictive model.

More specifically, in 240, for each value of f, one searches for the M-uplet of projectors w_fⁱ, i=1, . . . , M of respective dimensions I₁, I₂, . . . , I_M making it possible to minimise said penalised metric, as expressed in (6-1), (6-2), (6-3) in the case M=3. This minimisation takes place sequentially by ALS, projector by projector in 240₁, . . . , 240_M, the other projectors being considered as constant.

Advantageously, at each of these steps, the minimisation relates to the elements of the projector in question, as described in relation with the expressions (8-1) to (8-3) and (9-1) to (9-3).

More specifically, at each step, the elements of the projector making it possible to minimise the penalised metric are obtained by means of a thresholding of the elements of the projector minimising the non-penalised metric, as described in relation with the expressions (10-1), (10-2); (13-1) to (13-3); (16-1) to (16-3).

It is important to note that at each iteration f, only the penalisable elements are taken into account for the minimisation of the metric. To this end, one eliminates from the set of penalisable indices Ω_i,f, the indices k∈Ω_i,f for which the elements w_f,kⁱ have not been penalised at the iteration f, to generate the set of penalisable indices relative to the following iteration, Ω_i,f+1.

The tensor of the predictive model, B_u^f, is obtained from the projectors w_φⁱ, i=, . . . , M, Ω=1, . . . , f as in the conventional REW-NPLS regression method. Details of the calculation of B_u^f from the projectors in question will be found in appendix A of the article of A. Eliseyev et al. published in 2017, cited in the introduction part. The elements of the predictive model B_u^f corresponding to the indices for which the elements w_f,kⁱ have been penalised, are zeroed in 250.

After each iteration f, one eliminates in 260, from XX_u and XY_u the redundant information already taken into account in the regression on the latent variables space, at the preceding iterations 1, . . . , f. The cross-covariance tensor XY_u in which this redundancy has been eliminated is used for the new PARAFAC decomposition iteration. Step 260 is known by the name deflation in the NPLS regression.

At the end of the F iterations, in 270, one consequently has available F predictive models, B_u−1^f, β_u−1^f, f=1, . . . , F.

After the acquisition of new input calibration data, X_u+1 and output calibration data, Y_u+1 blocks at the following calibration step, 280, one determines in 290, the optimal predictive model, {B_u−1^f, β_u−1^f}. This predictive model is used to calculate the command tensor from the observation tensor, up to the following calibration step, u+2.

FIG. 3 schematically represents the mechanism for updating the predictive model in the calibration method of FIG. 2.

In 310 is represented the taking into account of a new input calibration data block, represented by the tensor X_u and a new output calibration block, represented by the tensor Y_u.

The covariance XX_u and cross covariance, XY_u, tensors are deduced therefrom in 320.

For each value of f, a PARAFAC decomposition of the cross-covariance tensor is carried out using an alternating least squares (ALS) method in 330, each step of the PARAFAC decomposition relating to a projector, this decomposition being regularised by means of a penalisation term relating to this projector. The regularisation results in a thresholding of certain elements of the projector obtained by the conventional least squares method, that is to say non-penalised.

At each PARAFAC iteration represented by the loop 335, one deduces therefrom a predictive model {B_u^f, β_u^f} in 340.

For each predictive model, one eliminates in 350 the redundant information in the cross-covariance tensor (deflation) and one carries out a new PARAFAC iteration on the covariance tensor XX_u and the new cross covariance tensor XY_u (noted XY_n+1 in the figure) deflated of said redundancy, 360.

The recursive validation loop is represented in 380. It makes it possible to determine among the predictive models {B_u−1^f, β_u−1^f}, f=1, . . . , F, obtained at the preceding calibration step, that which minimises the prediction error on the new input, X_u+1, and output, Y_u+1 calibration data blocks.

Claims

1. Method for calibrating on-line a direct neural interface intended to receive a plurality of electrophysiological signals during a plurality of observation windows associated with observation times, to form an input tensor and to deduce therefrom, by means of a predictive model, an output tensor providing command signals intended to command one or more effectors, said calibration is carried out during a plurality of calibration steps, calibration step u being carried out from input calibration data, represented by an input calibration tensor, Xu∈□(I1× . . . ×IN), and output calibration data represented by an output calibration tensor Yu∈□(J1× . . . ×JM), said calibration step implementing a REW-NPLS regression comprising: the calculation (220) of a covariance tensor XXu from a covariance tensor XXu−1 obtained during a preceding calibration step and the tensor Xu, as well as the calculation of a cross covariance tensor XYu from a cross covariance tensor XYu−1 obtained during said preceding calibration step and the tensors Xu and Yu;a PARAFAC iterative decomposition of the cross covariance tensor XYu, each PARAFAC iteration being associated with a dimension f of the latent variables space on which said regression is carried out, each iteration providing a plurality of projectors wfi, i=1, . . . , M of respective sizes I1, I2, . . . IM and a predictive model defined by a prediction coefficients tensor Buf∈□(I1× . . . ×IM)×(J1× . . . ×JM) and a bias tensor βuf∈□(J1× . . . ×JM);characterised in that each PARAFAC iteration comprises a sequence of M elementary steps (2401, 2401, . . . 240M) of minimisation of a metric according to the alternating least squares method, each elementary minimisation step relating to a projector and considering the others as constant, said metric comprising a penalisation term that is a function of the norm of this projector, the elements of this projector not being subjected to a penalisation during a PARAFAC iteration f not being penalisable during following PARAFAC iterations.

2. Method for calibrating on-line a direct neural interface according to claim 1, characterised in that an elementary minimisation step relating to a projector wfi aims to minimise the penalised metric

3. Method for calibrating on-line a direct neural interface according to claim 1, characterised in that an elementary minimisation step relating to a projector wfi aims to minimise the penalised metrics

4. Method for calibrating on-line a direct neural interface according to claim 3, characterised in that q=0 and that the elements wf,ki resulting from the minimisation are given by: (wf,ki)L0=0 if k∈Ωi,f and |(wf,ki)LS|≤ThL0i (wf,ki)L0=(wf,ki)LS if not

5. Method for calibrating on-line a direct neural interface according to claim 3, characterised in that q=½ and that the elements wf,ki resulting from the minimisation are given by: (wf,ki)L0=0 if k∈Ωi,f and |(wf,ki)LS|≤ThL0,5i (wf,ki)L0=arg min(h0,5(0),h0,5(Bi·(wf,ki)LS)) if k∈Ωi,f and |(wf,ki)LS|>ThL0,5i (wf,ki)L0=(wf,ki)LS if not

6. Method for calibrating on-line a direct neural interface according to claim 3, characterised in that q=1 and that the elements wf,ki resulting from the minimisation are given by: (wf,ki)L1=0 if k∈Ωi,f and |(wf,ki)LS|≤ThL1i (wf,ki)L1=sgn((wf,ki)LS)(|(wf,ki)LS|−ThL1i) if k∈Ωi,f and |(wf,ki)LS|>ThL1i (wf,ki)L0=(wf,ki)LS if not

7. Method for calibrating on-line a direct neural interface according to claim 1, characterised in that the electrophysiological signals are ECoG signals.