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 (240.sub.1, 240.sub.1, . . . 240.sub.M) 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.

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, X.sub.u∈□.sup.(I.sup.1.sup.× . . . ×I.sup.N.sup.), and output calibration data represented by an output calibration tensor Y.sub.u∈□.sup.(J.sup.1.sup.× . . . ×J.sup.M.sup.), said calibration step implementing a REW-NPLS regression comprising: the calculation (220) of a covariance tensor XX.sub.u from a covariance tensor XX.sub.u−1 obtained during a preceding calibration step and the tensor X.sub.u, as well as the calculation of a cross covariance tensor XY.sub.u from a cross covariance tensor XY.sub.u−1 obtained during said preceding calibration step and the tensors X.sub.u and Y.sub.u; a PARAFAC iterative decomposition of the cross covariance tensor XY.sub.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.sub.f.sup.i, i=1, . . . , M of respective sizes I.sub.1, I.sub.2, . . . I.sub.M and a predictive model defined by a prediction coefficients tensor B.sub.u.sup.f∈□.sup.(I.sup.1.sup.× . . . ×I.sup.M.sup.)×(J.sup.1.sup.× . . . ×J.sup.M.sup.) and a bias tensor β.sub.u.sup.f∈□.sup.(J.sup.1.sup.× . . . ×J.sup.M.sup.); characterised in that each PARAFAC iteration comprises a sequence of M elementary steps (240.sub.1, 240.sub.1, . . . 240.sub.M) 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 w.sub.f.sup.i aims to minimise the penalised metric $\underset{{\tilde{w}}_f^i}{\min }\left({.Math.{\underset{\_}{\mathrm{XY}}}_{u\left(i\right)}-{{\tilde{w}}_f^i\left(w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1\right)}^T.Math.}^2+{\lambda }_i{.Math.{\tilde{w}}_f^i.Math.}_{q,{\Omega }_{i,f}}\right)$ where {tilde over (w)}.sub.f.sup.i is a non-normalised projector and $w_f^m=\frac{{\tilde{w}}_f^m}{.Math.{\tilde{w}}_f^m.Math.},$ m=1, . . . M; m≠i are normalised projectors, XY.sub.u(i) is a matrix of size $I_i\times \underset{\underset{n\ne i}{n=1}}{\overset{N}{.Math.}}{I}_n$ obtained by deployment of the cross covariance tensor XY.sub.u while taking the mode i as reference, .Math. is the Kronecker product, λ.sub.i is a strictly positive coefficient, ∥w∥.sub.q,Ω.sub.i,f, q=0, ½, 1 is a norm defined by ${.Math.w_f^i.Math.}_{0,{\Omega }_{i,f}}=\underset{k\in {\Omega }_{i,f}}{.Math.}\left(1-{\delta }_{0,w_{f,k}^i}\right),{.Math.w_f^i.Math.}_{1/2,{\Omega }_{i,f}}=\underset{k\in {\Omega }_{i,f}}{.Math.}\sqrt{.Math.w_{f,k}^i.Math.},{.Math.w_f^i.Math.}_{1,{\Omega }_{i,f}}=\underset{k\in {\Omega }_{i,f}}{.Math.}.Math.w_{f,k}^i.Math.$ where w.sub.f,k.sup.i is the k.sup.th element of the projector w.sub.f.sup.i, δ is the Kronecker symbol, and Ω.sub.i,f is the set of indices corresponding to the penalisable elements of the projector w.sub.f.sup.i at the PARAFAC iteration f.

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 w.sub.f.sup.i aims to minimise the penalised metrics $\underset{w_{f,k}^i}{\min }\left({.Math.z_i^k-{w_{f,k}^i\left({\left(w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1\right)}^T\right)}^T.Math.}^2+{\lambda }_i{g}_q\left(w_{f,k}^i\right)\right),$ k=1, . . . , I.sub.i where z.sub.i.sup.k is the k.sup.th column vector of the matrix XY.sub.u(i) of size $I_i\times \underset{\underset{n\ne i}{n=1}}{\overset{N}{.Math.}}{I}_n$ obtained by deployment of the cross covariance tensor XY.sub.u while taking the mode i as reference, .Math. is the Kronecker product, λ.sub.i is a strictly positive coefficient, and the functions g.sub.q, q=0, ½, 1 are defined by:
g.sub.0(w.sub.f,k.sup.i)=1−δ.sub.0,w.sub.f,k.sub.i if k∈Ω.sub.i,f and g.sub.0(w.sub.f,k.sup.i)=0 if not
g.sub.1/2(w.sub.f,k.sup.i)=√{square root over (|w.sub.f,k.sup.i|)} if k∈Ω.sub.i,f and g.sub.1/2(w.sub.f,k.sup.i)=0 if not
g.sub.1(w.sub.f,k.sup.i)=|w.sub.f,k.sup.i| if k∈Ω.sub.i,f and g.sub.1(w.sub.f,k.sup.i)=0 if not where Ω.sub.i,f is the set of indices k corresponding to the penalisable elements w.sub.f,k.sup.i of the projector w.sub.f.sup.i at the PARAFAC iteration f.

4. Method for calibrating on-line a direct neural interface according to claim 3, characterised in that q=0 and that the elements w.sub.f,k.sup.i resulting from the minimisation are given by:
(w.sub.f,k.sup.i).sub.L0=0 if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|≤Th.sub.L0.sup.i
(w.sub.f,k.sup.i).sub.L0=(w.sub.f,k.sup.i).sub.LS if not where (w.sub.f,k.sup.i).sub.LS represents the minimisation solution in the least squares sense: ${\left(w_{f,k}^i\right)}_{\mathrm{LS}}=\frac{z_i^k\left(w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1\right)}{{.Math.w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1.Math.}^2}$ and that the threshold Th.sub.L0.sup.i is given by ${\mathrm{Th}}_{L0}^i=\frac{\sqrt{{\lambda }_i}}{.Math.w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1.Math.}.$

5. Method for calibrating on-line a direct neural interface according to claim 3, characterised in that q=½ and that the elements w.sub.f,k.sup.i resulting from the minimisation are given by:
(w.sub.f,k.sup.i).sub.L0=0 if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|≤Th.sub.L0,5.sup.i
(w.sub.f,k.sup.i).sub.L0=arg min(h.sub.0,5(0),h.sub.0,5(B.sup.i.Math.(w.sub.f,k.sup.i).sub.LS)) if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|>Th.sub.L0,5.sup.i
(w.sub.f,k.sup.i).sub.L0=(w.sub.f,k.sup.i).sub.LS if not where h.sub.0,5(w.sub.f,k.sup.i)=∥z.sub.1.sup.k−w.sub.f,k.sup.i(w.sub.f.sup.3≃w.sub.f.sup.2).sup.T∥.sup.2+λ.sub.i√{square root over (|w.sub.f,k.sup.i|)} and where (w.sub.f,k.sup.i).sub.LS represents the minimisation solution in the least squares sense: ${\left(w_{f,k}^i\right)}_{\mathrm{LS}}=\frac{z_i^k\left(w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1\right)}{{.Math.w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1.Math.}^2}$ and the threshold Th.sub.L0,5.sup.i is given by ${\mathrm{Th}}_{L0,5}^i=\frac{3}{4}{\left(\frac{{\lambda }_i}{{.Math.w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1.Math.}^2}\right)}^{\frac{2}{3}}.$

6. Method for calibrating on-line a direct neural interface according to claim 3, characterised in that q=1 and that the elements w.sub.f,k.sup.i resulting from the minimisation are given by:
(w.sub.f,k.sup.i).sub.L1=0 if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|≤Th.sub.L1.sup.i
(w.sub.f,k.sup.i).sub.L1=sgn((w.sub.f,k.sup.i).sub.LS)(|(w.sub.f,k.sup.i).sub.LS|−Th.sub.L1.sup.i) if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|>Th.sub.L1.sup.i
(w.sub.f,k.sup.i).sub.L0=(w.sub.f,k.sup.i).sub.LS if not where (w.sub.f,k.sup.i).sub.LS represents the minimisation solution in the least squares sense: ${\left(w_{f,k}^i\right)}_{\mathrm{LS}}=\frac{z_i^k\left(w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1\right)}{{.Math.w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1.Math.}^2}$ and the threshold Th.sub.L1.sup.i is given by ${\mathrm{Th}}_{L1}^i=\frac{{\lambda }_i}{.Math.w_f^M.Math.\mathrm{.Math.}{w}_f^{i+1}.Math.w_f^{i-1}.Math.\mathrm{.Math.}.Math.w_f^1.Math.}$

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

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0025] 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:

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

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

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

[0029] 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;

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

DETAILED DESCRIPTION OF PARTICULAR EMBODIMENTS

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

[0032] 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.

[0033] 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.

[0034] 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.

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

[0036] 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.

[0037] 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.sub.1× . . . ×I.sub.N.

[0038] 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.sub.1× . . . ×J.sub.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.

[0039] 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.

[0040] One will note hereafter X, the observation tensor at the time t. This tensor is consequently of order N and of dimension I.sub.1× . . . ×I.sub.N. It takes its values in a space X⊂□.sup.I.sup.1.sup.× . . . ×I.sup.N 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.sub.1× . . . ×J.sub.M. It takes its values in a space Y⊂□X.sub.u∈□.sup.J.sup.1.sup.× . . . ×J.sup.M.

[0041] 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.

[0042] If one notes {X.sub.u, Y.sub.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.sub.u−1 and XY.sub.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.sub.1, Y.sub.1}, {X.sub.2, Y.sub.2}, . . . {X.sub.u−1, Y.sub.u−1}.

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

XX.sub.u=γXX.sub.u−1+X.sub.u×.sub.1X.sub.u  (1-1)

XY.sub.u=γXY.sub.u−1+X.sub.u×.sub.1Y.sub.u  (1-2)

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

[0044] 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.

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

[0046] 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.sub.u. It will be recalled that, generally speaking, a PARAFAC decomposition of a tensor Z of order M, i.e. Z∈□.sup.I.sup.1.sup.×I.sup.2.sup.× . . . ×I.sup.M 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:

[00013]$\begin{array}{cc}\underset{\_}{Z}=\underset{r=1}{\overset{R}{.Math.}}{\theta }_r{w}_r^1\circ w_r^2\circ \mathrm{.Math.}\circ w_r^M+\underset{\_}{E}& \left(2\right)\end{array}$

with ∥w.sub.r.sup.m∥=1, m=1, . . . , M where ∘ is the external product and E is a residual tensor and .sub.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.

[0047] 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.sub.u a set of projection vectors, also designated projectors, w.sub.f.sup.1, . . . , w.sub.f.sup.M, of respective I.sub.1, I.sub.2, . . . , I.sub.M where the iteration number f gives the dimension of a latent variables space onto which the input tensor X.sub.u is projected.

[0048] 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.sub.u is equal to 3 (M=3), that is to say XY.sub.u∈□.sup.I.sup.1.sup.×I.sup.2.sup.×I.sup.3 with ∥XY.sub.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).

[0049] At the iteration f of the PARAFAC decomposition, one searches for the projection vectors w.sub.f.sup.1, w.sub.f.sup.2, w.sub.f.sup.3 which verify:

[00014]$\underset{}{\min }{.Math.{\underset{\_}{\mathrm{XY}}}_u-{\underset{\_}{Y}}_u.Math.}^2\mathrm{with}$${\underset{\_}{Y}}_u={\theta }_f{w}_f^1\circ w_f^2\circ w_f^3\mathrm{and}$$.Math.w_f^1.Math.=.Math.w_f^2.Math.=.Math.w_f^3.Math.=1$

[0050] 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.:

[00015]$\begin{array}{cc}\underset{w_f^1}{\min }{.Math.{\underset{\_}{\mathrm{XY}}}_{u\left(1\right)}-{w_f^1\left(w_f^3.Math.w_f^2\right)}^T.Math.}^2& \left(4\text{-}1\right)\\ \underset{w_f^2}{\min }{.Math.{\underset{\_}{\mathrm{XY}}}_{u\left(2\right)}-{w_f^2\left(w_f^3.Math.w_f^1\right)}^T.Math.}^2& \left(4\text{-}2\right)\\ \underset{w_f^3}{\min }{.Math.{\underset{\_}{\mathrm{XY}}}_{u\left(3\right)}-{w_f^3\left(w_f^2.Math.w_f^1\right)}^T.Math.}^2& \left(4\text{-}3\right)\end{array}$

[0051] where .Math. is the Kronecker product and Z.sub.(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.sub.(1) is defined by (z.sub.1.sup.1| . . . |z.sub.1.sup.I.sup.1)∈□.sup.I.sup.1.sup.×I.sup.2.sup.I.sup.3 where the Z.sub.1.sup.j are the columns of the matrix Z.sub.(1).

[0052] 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:

[00016]$\begin{array}{cc}{w}_f^1=\frac{{\underset{\_}{\mathrm{XY}}}_{u\left(1\right)}\left(w_f^3.Math.w_f^2\right)}{{.Math.w_f^3.Math.w_f^2.Math.}^2}& \left(5\text{-}1\right)\\ w_f^2=\frac{{\underset{\_}{\mathrm{XY}}}_{u\left(2\right)}\left(w_f^3.Math.w_f^1\right)}{{.Math.w_f^3.Math.w_f^1.Math.}^2}& \left(5\text{-}2\right)\\ w_f^3=\frac{{\underset{\_}{\mathrm{XY}}}_{u\left(3\right)}\left(w_f^2.Math.w_f^1\right)}{{.Math.w_f^2.Math.w_f^1.Math.}^2}& \left(5\text{-}3\right)\end{array}$

[0053] Each iteration f of the PARAFAC decomposition provides a plurality M of projectors w.sub.f.sup.1, . . . , w.sub.f.sup.M and a prediction coefficients tensor, or predictive model of the direct neural interface, defined by a prediction coefficients tensor B.sub.N.sup.f∈□.sup.(I.sup.1.sup.× . . . ×I.sup.N.sup.)×(J.sup.1.sup.× . . . ×J.sup.M.sup.) and a bias tensor β.sub.u.sup.f∈.sup.(J.sup.1.sup.× . . . ×J.sup.M.sup.), making it possible to predict the command tensor from the observation tensor. One thus obtains a plurality of models, B.sub.u.sup.f, f=1, . . . , F corresponding to different successive dimensions of the latent variables space.

[0054] One next determines, at the following calibration step, u+1, among B.sub.u.sup.f, f=1, . . . , F, the predictive model B.sub.u.sup.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.

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

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

[0057] At step 110, one estimates the vector w.sup.1.sub.f (corresponding to the epoch mode) from the matrix XY.sub.u(1) and the current vectors w.sub.f.sup.2 and w.sub.f.sup.3. This operation may be considered as a projection of the column vectors of the matrix XY.sub.u(1) onto the tensorial space generated by the modes 2 and 3.

[0058] At step 120, one estimates the vector w.sup.2.sub.f (corresponding to the frequency mode) from the matrix XY.sub.u(2) and the current vectors w.sub.f.sup.1 and w.sub.f.sup.3. This operation may be considered as a projection of the column vectors of the matrix XY.sub.u(2) onto the tensorial space generated by the modes 1 and 3.

[0059] At step 130, one estimates the vector w.sup.3.sub.f (corresponding to the spatial mode) from the matrix XY.sub.u(3) and the current vectors w.sub.f.sup.1 and w.sub.f.sup.2. This operation may be considered as a projection of the column vectors of the matrix XY.sub.u(3) onto the tensorial space generated by the modes 1 and 2.

[0060] Steps 110, 120, 130 are iterated until a convergence criterion is met, for example when the sum of the deviations between the vectors w.sub.f.sup.j of two successive ALS iterations is less than a predetermined threshold value.

[0061] 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.

[0062] 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:

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

where {tilde over (w)}.sub.f.sup.1, {tilde over (w)}.sub.f.sup.2 and {tilde over (w)}.sub.f.sup.3 are non-normalised projectors and

[00018]$w_f^1=\frac{{\tilde{w}}_f^1}{.Math.{\tilde{w}}_f^1.Math.},w_f^2=\frac{{\tilde{w}}_f^2}{.Math.{\tilde{w}}_f^2.Math.},w_f^3=\frac{{\tilde{w}}_f^3}{.Math.{\tilde{w}}_f^3.Math.}$

are the corresponding normalised projectors and λ.sub.1, λ.sub.2, λ.sub.3 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.sub.f.sup.i∥.sub.q,Ω.sub.i q=0, ½, 1 are defined by:

[00019]$\begin{array}{cc}{.Math.w_f^i.Math.}_{0,{\Omega }_{i,f},}=\underset{k\in {\Omega }_{i,f}}{.Math.}\left(1-{\delta }_{0,w_{f,k}^i}\right)& \left(7\text{-}1\right)\\ {.Math.w_f^i.Math.}_{1/2,{\Omega }_{i,f}}=\underset{k\in {\Omega }_{i,f}}{.Math.}\sqrt{.Math.w_{f,k}^i.Math.}& \left(7\text{-}2\right)\\ {.Math.w_f^i.Math.}_{1,{\Omega }_{i,f}}=\underset{k\in {\Omega }_{i,f}}{.Math.}.Math.w_{f,k}^i.Math.& \left(7\text{-}3\right)\end{array}$

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

[0063] 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 Ω.sub.i,j is higher.

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

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

where the vectors z.sub.1.sup.k, z.sub.2.sup.k, z.sub.3.sup.k are the k.sup.th column vectors of the matrices XY.sub.u(1), XY.sub.u(2), XY.sub.u(3) obtained by deployment of the cross covariance tensor XY.sub.u according to modes 1, 2 and 3, respectively.

[0065] The functions g.sub.q, q=0, ½, 1 apply to scalar values and are defined by:

g.sub.0(w.sub.f,k.sup.i)=1−δ.sub.0,w.sub.f,k.sub.i if k∈Ω.sub.i,f and g.sub.0(w.sub.f,k.sup.i)=0 if not  (9-1)

g.sub.1/2(w.sub.f,k.sup.i)=√{square root over (|w.sub.f,k.sup.i|)} if k∈Ω.sub.i,f and g.sub.1/2(w.sub.f,k.sup.i)=0 if not  (9-2)

g.sub.1(w.sub.f,k.sup.i)=|w.sub.f,k.sup.i| if k∈Ω.sub.i,f and g.sub.1(w.sub.f,k.sup.i)=0 if not  (9-3)

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

(w.sub.f,k.sup.i).sub.L0=0 if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|≤Th.sub.L0.sup.i  (10-1)

(w.sub.f,k.sup.i).sub.L0=(w.sub.f,k.sup.i).sub.LS if not  (10-2)

where (w.sub.f,k.sup.i).sub.LS represents the minimisation solution in the least squares sense (5-1), (5-2), (5-3), that is to say:

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

and the thresholds Th.sub.L0.sup.i, i=1, 2, 3 are given by:

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

[0067] 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.

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

(w.sub.f,k.sup.i).sub.L0=0 if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|≤Th.sub.L0,5.sup.i  (13-1)

(w.sub.f,k.sup.i).sub.L0=arg min(h.sub.0,5(0),h.sub.0,5(B.sup.i.Math.(w.sub.f,k.sup.i).sub.LS)) if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|>Th.sub.L0,5.sup.i   (13-2)

(w.sub.f,k.sup.i).sub.L0=(w.sub.f,k.sup.i).sub.LS if not  (13-3)

where h.sub.0,5(w.sub.f,k.sup.i)=∥z.sub.1.sup.k−w.sub.f,k.sup.i(w.sub.f.sup.3≃w.sub.f.sup.2).sup.T∥.sup.2+λ.sub.i√{square root over (|w.sub.f,k.sup.i|)} and the thresholds Th.sub.L0,5.sup.i, i=1, 2, 3 are given by:

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

[0069] The coefficients B.sup.i, i=1, 2, 3 are defined as the respective solutions of the cubics x(1−x).sup.2=C.sub.i, i=1, 2, 3, with:

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

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

(w.sub.f,k.sup.i).sub.L1=0 if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|≤Th.sub.L1.sup.i   (16-1)

(w.sub.f,k.sup.i).sub.L1=sgn((w.sub.f,k.sup.i).sub.LS)(|(w.sub.f,k.sup.i).sub.LS|−Th.sub.L1.sup.i) if k∈Ω.sub.i,f and |(w.sub.f,k.sup.i).sub.LS|>Th.sub.L1.sup.i   (16-2)

(w.sub.f,k.sup.i).sub.L0=(w.sub.f,k.sup.i).sub.LS if not  (16-3)

where the thresholds Th.sub.L1.sup.i, i=1, 2, 3 are given by:

[00025]$\begin{array}{cc}{\mathrm{Th}}_{L1}^1=\frac{{\lambda }_1}{.Math.w_f^3.Math.w_f^2.Math.}& \left(17\text{-}1\right)\\ {\mathrm{Th}}_{L1}^2=\frac{{\lambda }_2}{.Math.w_f^3.Math.w_f^1.Math.}& \left(17\text{-}2\right)\\ {\mathrm{Th}}_{L1}^3=\frac{{\lambda }_3}{.Math.w_f^2.Math.w_f^1.Math.}& \left(17\text{-}3\right)\end{array}$

[0071] 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.

[0072] 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.

[0073] At step 200, one takes into account new input and output calibration data, respectively in the form of an input tensor X.sub.u and an output tensor Y.sub.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.

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

[0075] At step 220, one updates the covariance XX.sub.u and cross covariance, XY.sub.u, tensors according to the relationships (1-1) and (1-2), while taking account of the forgetting factor.

[0076] 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.sub.u.sup.f∈□.sup.(I.sup.1.sup.× . . . ×I.sup.M.sup.)×(J.sup.1.sup.× . . . ×J.sup.M.sup.) and a bias tensor β.sub.u.sup.f∈□.sup.(f.sup.1.sup.× . . . ×J.sup.M.sup.).

[0077] At step 230, one initialises f with f=1 and the sets Ω.sub.i,1, i=1, . . . , M with Ω.sub.i,1={1, 2, . . . , I.sub.i} for i=1, . . . , M.

[0078] 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.

[0079] More specifically, in 240, for each value of f, one searches for the M-uplet of projectors w.sub.f.sup.i, i=1, . . . , M of respective dimensions I.sub.1, I.sub.2, . . . , I.sub.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.sub.1, . . . , 240.sub.M, the other projectors being considered as constant.

[0080] 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).

[0081] 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).

[0082] 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 Ω.sub.i,f, the indices k∈Ω.sub.i,f for which the elements w.sub.f,k.sup.i have not been penalised at the iteration f, to generate the set of penalisable indices relative to the following iteration, Ω.sub.i,f+1.

[0083] The tensor of the predictive model, B.sub.u.sup.f, is obtained from the projectors w.sub.φ.sup.i, i=, . . . , M, Ω=1, . . . , f as in the conventional REW-NPLS regression method. Details of the calculation of B.sub.u.sup.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.sub.u.sup.f corresponding to the indices for which the elements w.sub.f,k.sup.i have been penalised, are zeroed in 250.

[0084] After each iteration f, one eliminates in 260, from XX.sub.u and XY.sub.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.sub.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.

[0085] At the end of the F iterations, in 270, one consequently has available F predictive models, B.sub.u−1.sup.f, β.sub.u−1.sup.f, f=1, . . . , F.

[0086] After the acquisition of new input calibration data, X.sub.u+1 and output calibration data, Y.sub.u+1 blocks at the following calibration step, 280, one determines in 290, the optimal predictive model, {B.sub.u−1.sup.f, β.sub.u−1.sup.f}. This predictive model is used to calculate the command tensor from the observation tensor, up to the following calibration step, u+2.

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

[0088] In 310 is represented the taking into account of a new input calibration data block, represented by the tensor X.sub.u and a new output calibration block, represented by the tensor Y.sub.u.

[0089] The covariance XX.sub.u and cross covariance, XY.sub.u, tensors are deduced therefrom in 320.

[0090] 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.

[0091] At each PARAFAC iteration represented by the loop 335, one deduces therefrom a predictive model {B.sub.u.sup.f, β.sub.u.sup.f} in 340.

[0092] 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.sub.u and the new cross covariance tensor XY.sub.u (noted XY.sub.n+1 in the figure) deflated of said redundancy, 360.

[0093] The recursive validation loop is represented in 380. It makes it possible to determine among the predictive models {B.sub.u−1.sup.f, β.sub.u−1.sup.f}, f=1, . . . , F, obtained at the preceding calibration step, that which minimises the prediction error on the new input, X.sub.u+1, and output, Y.sub.u+1 calibration data blocks.