Articles
https://doi.org/10.1038/s42256-021-00306-1

Dimensionality reduction to maximize prediction
generalization capability
Takuya Isomura

1,2✉ and Taro Toyoizumi

1,3✉

Generalization of time series prediction remains an important open issue in machine learning; earlier methods have either
large generalization errors or local minima. Here, we develop an analytically solvable, unsupervised learning scheme thatextracts the most informative components for predicting future inputs, which we call predictive principal component analysis(PredPCA). Our scheme can effectively remove unpredictable noise and minimize test prediction error through convex optimization. Mathematical analyses demonstrate that, provided with sufficient training samples and sufficiently high-dimensional
observations, PredPCA can asymptotically identify hidden states, system parameters and dimensionalities of canonical nonlinear generative processes, with a global convergence guarantee. We demonstrate the performance of PredPCA using sequential
visual inputs comprising handwritten digits, rotating three-dimensional objects and natural scenes. It reliably estimates distinct hidden states and predicts future outcomes of previously unseen test input data, based exclusively on noisy observations.
The simple architecture and low computational cost of PredPCA are highly desirable for neuromorphic hardware.

P

rediction is essential for both biological organisms1–3 andmachine learning4–6. In particular, they both need to predictthe dynamics of newly encountered sensory inputs (that is, testdata) based on and only on knowledge learned from a limited number of past experiences (that is, training data). Generalization error
is a standard measure of the generalization capability of predicting
the future consequences of previously unseen input data, which isdefined as the difference between the training and test prediction
errors. It is thus crucial for organisms and machines to find a prediction strategy with a small generalization error, because otherwise
their predictions will fail because of overfitting to the training data.
Despite the importance of generalizing prediction, current
mainstream machine learning approaches have some limitations.
The approaches can be categorized into three major groups, andtheir limitations are summarized as follows: (1) The most basic
prediction strategy is to learn a direct mapping from past to future
inputs in the form of an autoregressive model (Fig. 1a). Althoughautoregressive models are simple to construct and guarantee global
convergence, their predictions contain a large generalization error
because the mapping from the observations to the prediction isoften redundant, leading to severe overfitting when the numberof training samples is limited7,8. Thus, to make accurate predictions, low-dimensional (that is, concise) representations should be
extracted from high-dimensional (that is, redundant) sensory data.
(2) A dimensionality reduction technique can be used to obtain aconcise representation;9 however, this is often achieved separatelyfrom the prediction step—for example, by first applying an autoencoder to reduce the dimensionality10,11 and then employing along short-term memory to predict the sequence12 (Fig. 1b). The
first autoencoding step—which provides a low-dimensional representation that minimizes the loss for reconstructing the current
input—is the most basic dimensionality reduction strategy. One
problem with this approach is that autoencoders may preferentiallyextract observation noise that is useless for prediction, owing to its
extra variance. From a prediction perspective, it is more helpful to

reduce the dimensionality to minimize the prediction error, similar to the approach used in time-lagged autoencoders (TAEs)13 andtheir variants14,15 (Fig. 1c). These approaches combine predictions
with dimensionality reduction in a single architecture. (3) A major
approach to time-series prediction is to construct a state-spacemodel (SSM). SSMs, which include the Kalman filter16 and its nonlinear variants17,18, simultaneously perform dimensionality reduction and prediction (Fig. 1d). From this model-based perspective,
the best prediction is achieved when an SSM employs the statesand parameters that match the true properties of the external system. However, the problem becomes difficult when both the hidden states and system parameters are unknown. In particular, theirpredictions become inaccurate owing to nonlinear interactions
between the uncertainties in hidden states and parameters, because
they can create spurious solutions. Furthermore, the dimensionalityof hidden states, which is essential for prediction accuracy, is difficult to optimize. Conventional model selection approaches using
some information criterion19–21, structural risk22 or cross-validation23
would fail to identify the optimal dimensionality when the state orparameter estimation converges to a suboptimal solution. In short,
all three approaches have essential drawbacks that interfere with thegeneralization of accurate predictions.
To overcome these limitations, we establish a method that cansolve this simultaneous optimization problem of hidden states, system parameters and dimensionality with a global convergence guarantee. We develop an unsupervised learning scheme for extracting
features that are essential for prediction, which we call predictive
principal component analysis (PredPCA). It is formally derived
from the minimization of the squared prediction error and can
extract low-dimensional predictive features from high-dimensional
sensory inputs, even in the presence of observation noise that ismuch larger than the signals themselves. This robustness is because
PredPCA conducts post hoc dimensionality reduction to extracta concise representation of the predicted input (Fig. 1e), unlike
autoencoders or SSMs. Moreover, the architecture of PredPCA isLaboratory for Neural Computation and Adaptation, RIKEN Center for Brain Science, Wako, Japan. 2Brain Intelligence Theory Unit, RIKEN Center for
Brain Science, Wako, Japan. 3Department of Mathematical Informatics, Graduate School of Information Science and Technology, The University of Tokyo,
Tokyo, Japan. ✉e-mail: takuya.isomura@riken.jp; taro.toyoizumi@riken.jp
1434

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCEab

Autoregressive modelϕt
st+1∣t

st

ut

c
s~t

~
~
st+1∣t – εt+1∣tut

ut+1∣tResults

Overview of PredPCA. In this work, we assume that hidden states
xt generate higher-dimensional sensory inputs st as follows:
st = g (xt ) + ωt ,

d

TAE
~
st+1∣tdimensions and number of past observations used for prediction,and can provide the global minimum of the test prediction error.
We mathematically and numerically demonstrate that filtering out
unpredictable noise by using PredPCA is essential to maximizing
the prediction generalization capability.

Autoencoder (PCA)s t – εtSSMst+1∣tstand the dynamics of hidden states are described by
xt+1 = f (xt , xt−1 , xt−2 , …) + zt ,u~t+1∣text∣txt+1∣tPredPCAϕt

st+1∣t – εt+1∣t

st+1∣tut+1∣tFig. 1 | Five different prediction model structures. A black bar denotes
a layer of a neural network, whereas blue and green trapezoids denotesynaptic weight matrices for prediction and dimensionality reduction,
respectively. a, Naive autoregressive models directly compute themaximum likelihood estimator of the next input st+1|t based on thebases ϕt ≡(sTt , sTt−1 , …, sTt−Kp +1)Tthat summarize current and pastobservations. b, Two-step prediction models first extract a concise
representation ut using an autoencoder (or principal component analysis,
PCA) by minimizing the loss εt and then predict the next representation
ut+1|t using a recurrent neural network. c, TAEs and their variants combine
predictions with dimensionality reduction by performing eigenvalue
decomposition or singular value decomposition of the transition matrixof the input. Note that s̃t ≡ Σs−1/2 st denotes the normalized inputs,s̃t+1|t ≡ Σs−1/2 st+1|t denotes the normalized predicted inputs, and Σsdenotes the actual input covariance. d, SSMs update the hidden state
estimator xt|t based on the previous state and current input, and predict
the next state xt+1|t and input st+1|t. e, PredPCA first computes the maximum
likelihood estimator st+1|t based on multi-timestep basis functions ϕt andthen extracts a concise representation ut+1|t, by minimizing the prediction
error εt+1|t. This scheme can effectively filter out the causes of thegeneralization error.

suitable for noise reduction because it predicts the subsequent
input based on multi-timestep basis functions, unlike TAEs andtheir variants. These properties allow PredPCA to find hidden
states (refer to blind source separation24) and to perform long-term
prediction reliably and accurately. In particular, system parameter
identification25,26 using PredPCA contrasts with conventional methods. It is guaranteed to asymptotically identify the true parametersof canonical nonlinear systems (see below for the definition) in thelarge sample-size limit, when the mappings from hidden states tosensory inputs are sufficiently high-dimensional. In addition, basedon Akaike’s statistics19,27, we analytically derive a mathematical formula that estimates the test prediction error of PredPCA. It showsthat the generalization error is proportional to an entropy that is dueto the sampling fluctuation27. The minimization of this formula can
optimize unknown free parameters, including the rank of system

(1)

(2)

where zt and ωt are mutually independent white noises, with zeromeans and covariances Σz and Σω (Fig. 2a, left, and Methods section ‘System'). Process noise zt adds stochasticity into the hidden
state dynamics, while observation noise ωt represents any unpredictable fluctuations that we wish to remove. PredPCA is applicable
to sensory data generated from systems involving either Gaussianor non-Gaussian noise to extract features and characterize systemproperties, although the identification of noise distributions is morestraightforward when provided with Gaussian noise (see Methods
section ‘System parameter identification'). Although this paper
focuses on white noise, PredPCA’s outcomes would also be accuratewith coloured noise as long as the auto-correlation time constant of
ωt is smaller than that of xt. To apply PredPCA to continuous-timesystems, the time bin size should be determined depending on theirapplications. Table 1 presents the glossary of expressions.
PredPCA aims to extract the components containing the mostinformation for predicting the next input st+1 based on current andpast observations st , st−1 , …, st−Kp +1. With this in mind, we consider a linear neural network whose output is given by
ut+1|t = Vϕt ,(3)where ut+1|t is an Nu-dimensional vector of encoders, V is a (horizontally long) Nu × Nϕ encoding synaptic weight matrix, and()T
ϕt ≡ sTt , sTt−1 , …, sTt−Kp +1 is an Nϕ-dimensional vector of linear

basis functions that summarize current and past observations. We
refer to this linear encoder as a neural network, intending an analogy
to biological neural networks and to highlight potential applications
to neuromorphic computation (see ‘Discussion' section for further
details). Unlike standard PCA28,29 and autoencoders10,11, which minimize the reconstruction error in the current input, PredPCA minimizes the prediction error ϵt+1|t ≡ st+1 − WT ut+1|t , defined as thedifference between the actual next input at t + 1 and the prediction
based on inputs up to t. Here, WT is an Ns × Nu decoding synaptic
weight matrix used for predicting the next input st+1 based on the
concise encoders ut+1|t (where we introduced the transposed matrix
WT rather than W for a notational reason that will become clear
below). PredPCA’s cost function L is defined as the expectation of
the squared prediction error over the training period T:
L≡

2 
1 
ϵt+1|t 
2
q

(4)

1 ∑T
Here, ⟨•⟩q ≡ T t=1 • indicates the expectation over the empirical distribution q. By minimizing this cost function with respect
to V, we obtain the optimal encoding weights as V = WQ, where
⟨
⟩ ⟨
⟩− 1
Q ≡ st+1 ϕTt q ϕt ϕTt q (see Methods section ‘Derivation of

PredPCA'). Thus, ut+1|t = Wst+1|t holds, where st+1|t = Qϕt is the

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

435

Articles

NATurE MACHInE InTEllIGEnCE
t

t+1

x

x

x

g
Sensory input

s

s

s

0

u

u

–6
–6

0

u

u

xt+1∣t,10

0

HMM
0

50

6

0

0

6

–6

0

Error (%)

16

d

100

12

60

0
103 104 105
Samples

8

A
B

Σx
Σω
Σz

4
0
103

100

104

Categorization error (%)

6

xt+1∣t,9

Test prediction error (%)

PCA of ϕt+ICA
TICA+ICATAE+ICADMD+ICADMD of ϕt+ICA
SSM(KF/BF)+ICA

0
xt+1∣t,5–6–6

cParameter estimation error (%)

0.0012 %1.2 %PredPCA+ICA
PCA+ICA
Autoencoder+ICA–6

xt+1∣t,3

xt+1∣t,7

Ascending
Fibonacci

6

6

–6

b

0

s

W
Encoders

66s

xt+1∣t,8s

0
xt+1∣t,1

Q
s

0

–6

–6

s

–6

Predicted input

6x
xt+1∣t,2

Hidden states

6

6

f

xt+1∣t,6

t–1

105

Error (%)

t–2

xt+1∣t,4

a

51
50
49
10

55

Nu

40

50

103

Training samples

104

105

Training samplese
Ascending
Fibonacci
1

2

3

4

5

6

7

8

9

10

99,99299,994

99,99699,998105

Prediction steps

Fig. 2 | PredPCA of handwritten digit sequences. a, On the left we show the system, comprising a generative process (top) and a neural network that
follows PredPCA (bottom, shaded). The network is trained with an image sequence st of handwritten digits generated from the dynamics of 10-dimensional
hidden states xt, each element of which expresses one of the ten digits. On the right, we show the 10-dimensional independent encoders (that is, hidden
state estimator) xt+1|t obtained using PredPCA and ICA. 2 × 104 test samples that are colour-coded by their digit are plotted. b, Comparison with related
methods in terms of the mean categorization error (that is, false discovery rate), obtained by averaging categorization errors of ten elements of xt+1|t. The
digits are introduced in ascending order (blue) and Fibonacci sequence (red). An SSM based on a Kalman filter (KF) is used for the ascending sequence,
while that based on a Bayesian filter (BF) is used for the Fibonacci sequence. The green bars indicate the minimum categorization error among 20 different
realizations of digit sequences. c, Parameter estimation error measured by the squared Frobenius
norm() ratio, where the difference between the ground truth
parameter matrix θ and its estimator θ is divided by their norm, error = |θ − θ|2F / max |θ|2F , |θ|2F . We assume here that the ascending-order handwritten
digit sequence is generated from a linear system comprising st = Axt + ωt and xt+1 = Bxt + zt (A, black, B, red). The covariance matrices Σx (blue), Σω (green),and Σz (grey, inset) are associated with xt, ωt and zt, respectively. d, Test error in predicting
in the ascending-order2  digit images
 the next handwrittensequence, measured by the normalized mean squared error over test samples, error = st+1 − WT ut+1|t  / |st+1 |2 . The red line (in the main andinset panels) represents a theoretical prediction obtained using equation (7). The blue line denotes the lower bound of the error, calculated via supervised
learning. The inset depicts the dependence of the test prediction error on the encoding dimensionality Nu (when T = 6,000), where Nu = 10 (green line) is
optimal. b−d are obtained with 20 different realizations of digit sequences and the error bars indicate the standard deviation, although some error bars ind are hidden by the circles. e, Long-term prediction using PredPCA and ICA. A winner-takes-all operation is applied to make greedy predictions of the digit
sequences. After receiving the first 40 digits, unless those initial digit images are outliers, the network can predict the next 105 digits (and more) without
any categorization error. See Supplementary Methods 1 and 2 for further details.

maximum likelihood estimator of st+1. The synaptic weight matrix
W is updated by gradient descent on L. After some additionaltransformations (see Methods section ‘Derivation of PredPCA'),
we obtain
⟨
(
)T ⟩
∂L
T
Ẇ ∝ −= ut+1|t st+1 − W ut+1|t
(5)
∂W
q
436

The fixed point of equation (5) yields the transpose of optimal
decoding weights that minimize L. The solution ensures that theencoders ut+1|t achieve the optimal representation for prediction.
Equation (5) is equivalent to the subspace rule of PCA28, exceptthat ut+1|t encodes the future state at time t + 1 instead of the state
at time t (that is, the standard PCA uses ut|t). This means that
PredPCA, which is defined by the prediction error minimization,

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCE
Table 1 | Glossary of expressions
Expression

Descriptionst

Observation

ψt

Hidden bases

xt

Hidden states

ωt

Observation noisezt

Process noiseA

Observation matrix (gt = Aψt)

B

State transition matrix (ft = Bψt)Σs, Σψ, Σx, Σω, Σz

Covariance matrices of st, ψt, xt, ωt, zt

NsDimensionality of observation

Nψ

Dimensionality of hidden bases

Nx

Dimensionality of hidden states

ut+k|tEncoders

ϕtBasis functions

V

Encoding synaptic weight matrix

W

Transpose of decoding synaptic weightmatrixNu

Dimensionality of encoders

NϕDimensionality of basis functions〈∙〉q

Expectation over empirical distribution q〈∙〉

Expectation over true distribution p

can be decomposed into two steps: computing the maximum likelihood estimator of st+1, st+1|t, followed by a post hoc PCA of st+1|t using
the eigenvalue decomposition (Fig. 1e). Owing to the global convergence property of the subspace rule for PCA30, the global convergence of equation (5) is also guaranteed. In essence, PredPCA is aconvex optimization. Crucially, however, only PredPCA (but not thestandard PCA) can effectively filter out unpredictable observationnoise, as we demonstrate numerically below and mathematically in
Methods section ‘Filtering out observation noise'. It is straightforward to extend PredPCA to multi-step predictions (see Methods
section ‘Derivation of PredPCA' for further details). We note that
although this paper focuses on the prediction of subsequent inputs(that is, autoregression), it is straightforward to apply PredPCA tominimize the generalization error for a class of regression tasks. The
formulation for this is performed simply by supposing that the hidden states xt generate both observations st and a high-dimensionaltarget signal yt and by replacing the prediction error εt+1|twith εt ≡ yt − WTVϕt.
After extracting the hidden states by using PredPCA, weemploy independent component analysis (ICA)31,32, which canseparate the extracted states into independent components aslong as the true hidden states of the external milieu are actuallymutually independent. For example, when the network observes asequence of handwritten digits generated using the MNIST dataset33, PredPCA followed by ICA generates 10-dimensional independent encoders xt+1|t, each element of which encodes one of
the ten possible digits (Fig. 2a, right). The detailed procedure to
extract xt+1|t from ut+1|t is provided in Methods section ‘Asymptotic
linearization theorem'.
Previous works have developed methods combining futuredata predictions with dimensionality reduction, for example,
time-lagged independent component analysis (TICA)14, TAE13
and dynamic mode decomposition (DMD)15. When ϕt = st,
PredPCA is involved in this family of methods—thus, one
may view PredPCA as a combination of these methods andautoregressive models based on high-dimensional, multi-timestep
basis functions. This construction enables PredPCA to effectively filter out observation noise and reduce test prediction error(see below).
Key analytical discoveries. We conducted comprehensive mathematical analyses to rigorously demonstrate the performance andstatistical properties of PredPCA. In particular, we demonstrated
the following two key properties. First, it is mathematically guaranteed that PredPCA can identify the optimal (explained below)
hidden state representation and parameter estimators—up to alinear transformation that does not affect prediction accuracy—
for general linear systems and, asymptotically, even for nonlinear
systems (Methods sections ‘Asymptotic linearization theorem' and‘System parameter identification'). When equations (1) and (2) are
involved in a class of canonical nonlinear systems defined by equations (8) and (9), a set of hidden states, parameters and dimensionalities that characterize a system is uniquely determined up toa trivial linear ambiguity (Methods section ‘System'). Under thiscondition, while using a linear neural network for the encoding,the asymptotic linearization theorem34 ensures that PredPCA will
extract the true hidden states when the hidden state dimensionality is large and the input dimensionality is sufficiently larger thanthe hidden state dimensionality. Briefly, this is because projecting the high-dimensional input onto the directions of the major
eigenvectors of the input covariance effectively magnifies the linearly transformed components of the hidden states included in the
input, while filtering out the nonlinear components (see Methods
section ‘Asymptotic linearization theorem' for its mathematicalstatement and the conditions for application; see ref. 34 for themathematical proof).
Owing to this linearization property, the hidden state estimator xt+1|t obtained using PredPCA asymptotically converges to alinear transformation of the maximum likelihood estimator of
⟨
⟩ ⟨
⟩− 1
hidden states xt+1, that is, xt+1 ϕTt q ϕt ϕTt q ϕt . Hence, PredPCA
provides the optimal hidden state representation for prediction.
Furthermore, the analytical expressions of the system parameter
estimators are derived as functions of xt+1|t, with a convergenceguarantee to the true parameter values in the large sample-sizeand system-size limits. These parameter estimators are calculated
by a simple iteration-free computation summarized in Table 2 andMethods section ‘System parameter identification'. In essence,
provided with sufficient but finite training samples, PredPCA canidentify the hidden states and parameters of large-scale canonical systems up to a small estimation error. This result is surprising because the reliable identification of the optimal hidden states
and the true parameters were previously only described within
the framework of supervised learning, whereas PredPCA can provide them by unsupervised learning without relying on the true
hidden states xt.
Second, PredPCA can maximize the prediction generalization
capability by minimizing the test prediction error
Ltest ≡

2 1 
ϵt+1|t 2(6)Here, 〈∙〉 ≡ ∫∙p(ϕt,st+1)dϕtdst+1 indicates the expectation over thetrue distribution p(ϕt,st+1) (note the difference from equation (4)).
In practice, however, the true distribution is unknown for a learner.
Thus, one needs to estimate equation (6) based on and only on
parameters estimated from the training data. In the framework of
the maximum likelihood estimation or squared error minimization,
the expectation of the test error is expressed as an Akaike information criterion (AIC)19 or network information criterion (NIC)20,
respectively. Similar to the derivation of AIC and NIC, we explicitly

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

437

Articles

NATurE MACHInE InTEllIGEnCE

Table 2 | Definitions and analytical solutions of estimators
Estimator

Definitionst+k|t⟨Analytical solution⟩ −1
ϕt ϕTt q ϕt

ψ t + k| t

st+k ϕTt q
PTs st+k|t

xt+k|t1/2 TΛ−P ψ ψ t + k| t
ψ

A

⟩ ⟨

⟩
(
)
st+k ϕTt Σϕ−1 ϕt + O T−1/2
⟨
⟩
(
)
−
1
Ωψ ψ t+k ϕTt Σϕ ϕt + O T−1/2⟨⟩()Ωx xt+k ϕTt Σϕ−1 ϕt + O T−1/2 + O (σ x )()
AΩψ−1 + O T−1/2()Ωx BΩψ−1 + O T−1/2 + O (σ x )

⟨

P

Ψ

⟨s
⟩ ⟨
⟩ −1
ψ t+1|t ψ Tt+1|t
xt+2|t ψ Tt+1|t
q
q
⟨
⟩ ⟨
⟩ −1
ψ t+1|t+1 ψ Tt|t
ψ t+2|t+2 ψ Tt|t

Σs

⟨

B

q

Ωψ Ψ Ωψ−1 + O T−1/2 + O (σ ψ )(qΣs + O T−1/2

)

Σψ

⟩
st sTt q
(
)
⟨
⟨
⟩
⟩
1
Ψ−1 ψ t+1|t+1 ψ Tt|t + ψ t|t ψ Tt+1|t+1 Ψ−T2Ωψ Σψ ΩψT + O T−1/2 + O (σ ψ )ΣxIΣx ≡ IΣωΣs − AΣψAΣω + O T−1/2 + O (σ ψ )ΣzΣx − BΣψBT()Ωx Σz ΩxT + O T−1/2 + O (σ x )Nψargmin L

Converge to Nψ when T > Tcψ

Nx

Largest spectrum gap of ΣψConverge to Nx when T > Tcx and σ x < σ cx

q

q

T

(

(

(

Nu

)
)

)

The external system is characterized by st = Aψt + ωt and xt+1 = Bψt + zt. Throughout the Article, a bold variable (for example, st+k|t) indicates the estimator of the corresponding italic variable (for example,
(
)
(√
)
1/2
and Σψ, respectively. The full-rank square matrix Ωψ and orthogonal matrix Ωx are ambiguity factors. O (σ x ) = O
andst+k). Ps and Pψ are sets of the major eigenvectors of ΣPred
Nx /Nψ + O N−
s
x
(
)
1/2
are linearization errors, where σx = σψ = 0 for any linear system. Tcx , Tcψ < ∞ are finite large constants and σ cx > 0 is a small positive constant. Refer to Methods for further details.
O ( σ ψ ) = O N−ψcompute the expectation of equation (6), with the optimized synaptic weights, asL≡ E{q} [Ltest ]Test error expectation

=



1
Ps
tr [Σs ] − tr PTs ΣPred
s
2


Training error

+

(7)

 
 3
Nϕ  T 
Ps + O T− 2
tr Ps Σs − ΣPred
s
2T


Generalization error

The derivation is presented in Methods section ‘Test prediction
error minimization'. Here, T is the number of training samples,
Ps is the first-to-Nu-th major eigenvectors of the predicted input⟨⟩≡ st+1|t sTt+1|t q (where WT W = Ps PTs holds atcovariance ΣPreds⟨ T⟩
the fixed point of equation (5)), and Σs ≡ st st q is the actual input
covariance. The expectation E{q}[∙] is taken over different empirical
distributions q, each of which comprises T training samples and isused to optimize synaptic weights.
The expectation of the test prediction error L is characterizedby two free parameters: the rank of encoding dimensions (Nu) andthe number of past observations used for the maximum likelihood
estimation (Kp), where Nϕ = KpNs. Optimizations of Nu and Nϕ performed while updating the synaptic weight using equation (5) provide the global minimum of L. The optimal encoding dimensionality
is guaranteed to converge to the true hidden basis dimensionality of
the canonical system for a large but finite T (Methods section ‘Test
prediction error minimization'). The second term of L, referred toas the generalization error, is associated with an entropy that is dueto the sampling fluctuation27. This term indicates that only the prediction error projected to the major eigenspace causes the generalization error, which highlights the importance of dimensionality
reduction to reduce the test prediction error. In short, naively minimizing the training error by using a large encoding dimensionality,
such as in autoregressive models, leads to overfitting; in contrast,
438

minimizing L provides the best encoding dimensionality and number of past observations to generalize the prediction.
For further details, please see Methods and Supplementary
Information. The aforementioned analytical results are empirically validated though numerical simulations by confirming thereliable identification of properties of canonical systems defined in
Methods section ‘System' (Supplementary Fig. 1a–c). Furthermore,
empirical observations imply that the outcomes of PredPCA can beutilized to identify the properties of more general classes of systems(for example, a class involving a Lorenz attractor; Supplementary
Fig. 1d–f), although system parameter identification beyond theclass defined in Methods section ‘System' has not yet been proved
mathematically. In what follows, we demonstrate the performanceof PredPCA using sequential visual inputs comprising handwritten digits, rotating three-dimensional (3D) objects and natural scenes (refer to Supplementary Methods sections 1 and 2 forsimulation protocols).
PredPCA provides optimal representation and parameters forprediction. In the first experiment (Fig. 2), we trained a neural network with MNIST handwritten digit images33 in ascending orderand in the Fibonacci sequence, wherein only the last digit was presented; however, these sequences involve some additional stochasticity (which corresponds to process noise zt) such that a digit was
replaced by a random one and a monochrome inversion occurredwith a small probability at each step (analogous to large noise thatinterferes with weak signal measurements: for example, movement
artefacts in electroencephalogram recordings35). In both cases,
PredPCA successfully extracted 10-dimensional features underlying the image sequences as they were relevant to predicting thesequences. The following ICA32 separated the extracted components
into independent hidden states. Each of the ensuing encoder neurons (that is, independent components, xt+1|t) selectively responded
to one of the ten digits without being taught their labels, as we
can see for the encoders trained with the ascending sequence in
Fig. 2a (right).
Irrespective of the sequence types (ascending order andFibonacci sequence), PredPCA and ICA precisely separated the digits into ten clusters in 10 dimensions with an average categorization

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCE
error of less than 2% (scored by false discovery rate; Fig. 2b). During
this process, PredPCA ignored any within-class differences in thedigit images that do not predict the next image (which correspondto observation noise ωt). Hence, PredPCA’s policy of dimensionality reduction to minimize the prediction error distinguishes it fromstandard PCA28,29 and autoencoders10,11—because PCA and autoencoders minimize the reconstruction error for the current input
st and thus preferentially extract the within-class differences in thedigit images owing to their extra variances. Even when the standard
PCA was applied to the past-to-current input sequence (that is, ϕt),
it failed to separate the digits because the hidden representation of
ϕt included more than 10-dimensional state space and thus the firstten major components of ϕt did not match the true hidden states
xt. Although the categorization errors of TICA14, TAE13 and DMD15
were smaller than those of PCA and an autoencoder, they still failed
to categorize some digits. This is because the former methods use
only a single step (that is, ϕt = st) to predict subsequent digit images(st+1). Moreover, even when using (st, st-1) or ϕt to predict (st+1, st) or
ϕt+1, they failed to categorize digits, because the extracted features
do not match the true hidden states (these results are similar to the
PCA of ϕt). The performance of the SSM and hidden Markov modelwith 10-dimensional state spaces was also poor because their larger
parameter estimation errors led them to a spurious solution or
local minimum.
In addition to accurate source separation, PredPCA could provide the optimal system parameters for the prediction (Fig. 2c).
These parameter estimators were computed simply by following
the definitions in Table 2. The differences between the parameter estimators obtained by PredPCA and those obtained bysupervised learning converged to zero as the number of training
samples increased, as predicted theoretically (Methods section‘System parameter identification'). These results corroborated that
PredPCA-based system parameter identification was applicable tosystems involving non-Gaussian noise. Consequently, the outcomesof PredPCA could reliably identify the transition rules underlying the ascending order (Extended Data Fig. 1a) and Fibonaccisequences (Extended Data Fig. 1b) in an unsupervised manner.
In essence, we demonstrated that each encoder obtained using
PredPCA corresponds to a digit, and the obtained state transition
matrix represents the estimated dynamics of digit sequences, which
can assign the meaning to these model parameters. These results
indicate the interpretability of PredPCA as the obtained model can
provide an explanation of the manner that the hidden dynamics
generate the sensory input.
The above outcomes allowed PredPCA to predict subsequent
digits reliably and accurately (Fig. 2d). Here, we see that although
PredPCA did not observe the hidden states directly, its test prediction error converged globally—with increasing training samples—tothe lower bound of the test prediction error computed via supervised learning that explicitly used the true hidden states for training.
This is as theoretically predicted by equation (7). Moreover, equation (7) successfully identified the optimal encoding dimensionality that minimized the test prediction error as Nu = 10, which also
matched the true hidden state dimensionality (Fig. 2d, inset). Thesematchings hold even in the absence of random replacement and/or monochrome inversion of digit images (Extended Data Fig. 1c).
Numerical observations indicate that PredPCA can reduce errors in
categorization, system identification, and prediction as the numberof past observations used for prediction (Kp) increases until reaching its finite optimum (Extended Data Fig. 1d). In contrast, linear
TAE and SSM (same as PredPCA with ϕt = st) failed to identify
the system properties, and thus generated a larger prediction error(Extended Data Fig. 2).
In particular, the long-term prediction of subsequent digits
highlights the virtue of PredPCA’s categorization and system identification accuracy—provided with a winner-takes-all operation,

the outcomes of PredPCA could recursively predict the subsequent digits without categorization errors for more than 105 steps(Fig. 2e). These results were minimally influenced by the assumed
transition mapping structures and training history (Extended Data
Fig. 3a–c), and the optimal model structure could be determinedthrough model selection based on the standard AIC (Extended
Data Fig. 3d). In contrast, SSMs tended to fail the long-term prediction depending on initial conditions and training history, eventhough they were provided with the winner-takes-all operation(Extended Data Fig. 3e).
PredPCA filters out observation noise and minimizes test prediction error. Next, the noise reduction and prediction generalization capabilities of PredPCA were examined using natural videos.
We trained a neural network by using images of 3D objects rotating
anti-clockwise36 as the input (Fig. 3a, furthest left). In short, the task
was to predict the opposite side of test object images (200 objects)
by observing only a half side of the images, based on the transition (that is, rotational) mapping learned from different training
object images (up to 800 objects). Here, we used the optimal linear bases ϕt to maximize PredPCA’s generalization capability (seeSupplementary Methods section 3 for the procedure). The abilityof PredPCA was experimentally confirmed by its successful predictions of the 30°–150° rotated images of previously unseen test
objects (Fig. 3a, middle row; see Supplementary Video 1 for predictions of 90° rotated images, where the right-hand-side images
are the predictions of the corresponding ground truth images on
the left-hand side).
In general, features extracted using PredPCA comprise categorical features that represent what the input is as well as dynamical features that express how it is moving. As the asymptotic linearization
theorem implies that the obtained encoders xt+k|t are linear superpositions of hidden states, (we define categorical features
as the aver)age of estimators, x̄t ≡ xt+30|t + … + xt+150|t /5, and dynamical
features as the deviation from their average, Δxt+k|t ≡ xt+k|t − x̄t .
Applying ICA to x̄t separated the categorical features into a sparse
representation, each dimension of which expresses a feature ofobjects (Fig. 3b, top). Applying an additional PCA to the dynamical features (for example, Δxt+90|t) provided the angle of 3D objects
as the first principal component (PC1; Fig. 3c, left). Although thecoordinate of the attractor changes depending on its category,
this treatment makes it easier to interpret the dynamics of hidden
states. We also observed the same property for xt+30|t, xt+60|t, xt+120|t
and xt+150|t.
Notably, these prediction and feature extraction capabilities
were largely retained even in the presence of an artificially added
large (white Gaussian) observation noise whose variance had thesame magnitude as the variance of original images, demonstrating the robustness of PredPCA’s outcomes (Fig. 3a,b, bottom, andFig. 3c, right; see Supplementary Video 2 for predictions of 90°
rotated images). The sampling fluctuation caused by the observation noise disturbed the prediction of minor components, and thuschanged the optimal encoding dimensionality (Fig. 3d).
We confirmed an earlier decrease of the test prediction error for
PredPCA relative to the naive autoregressive model as the numberof training samples increases (Fig. 3e). PredPCA generated a smaller
test prediction error relative to TICA, TAE, DMD and SSMs basedon the Kalman or Bayesian filters (Fig. 3f). These results indicatethat PredPCA could determine a plausible rule for rotating genericobjects. Remarkably, owing to the convex optimization, features
extracted using PredPCA are uniquely determined for any given
training dataset (even if the true system is unknown). This contrasts
with TAE and SSM, because their extracted features change depending on the initial conditions, order of supplying mini batches, orlevel of observation noise, even though they are trained with thesame dataset (Extended Data Fig. 4).

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

439

Articles

NATurE MACHInE InTEllIGEnCEabGround truth

Hidden state estimator xt

Predictions (without large noise)

Latest input

xt trained with large noise
Predictions (with large noise)

c

+60°

PC1 of ∆xt+90∣t

1
0
–1

+90°

+120°

With large noise

1

Activity

Activity

+30°

0
–1

0

90

180 270 360

Angle (°)

090

180 270 360

Angle (°)

+150°

d

e

Without noise

200

Optimal Nu

0°

Test prediction error (%)

Latest input

With noise

100

0
200

400

600

Training samples

800f
PredPCA
Autoregressive
model

120
80
40
0
200

400

600

Training samples

800

PredPCA
Autoregressive
model
TICA
TAE
DMD
SSM(KF)
SSM(BF)
0

20

40

60

Test prediction error (%)

Fig. 3 | PredPCA-based de-noising, hidden state extraction and subsequent input prediction of videos of rotating 3D objects. a, Snapshots of theprediction results. Latest input image (furthest left) and ground truth (top) and predicted images after 30°, 60°, 90°, 120° and 150° rotations, without(middle row) and with (bottom) artificially added observation noise. b, Images corresponding to 20-dimensional sparse representations ( x̄t) each
expressing a categorical feature of objects. These images were obtained (by applying ICA with super-Gaussianprior distribution to the first 20 principal)
components of PredPCA, averaged over different prediction points x̄t = xt+30|t + … + xt+150|t /5 (see Methods section ‘Asymptotic linearization
theorem' for the detail). These images visualize linear mappings from each independent component to the observation. c, Rotation of objects encoded in
the first principal component (PC1) of the dynamical features of xt+90|t (that is, Δxt+90|t = xt+90|t − x̄t). Here, the neural activity predicted the angle of90°-rotated future images, indicating that when observing an asymmetric object (as opposed to a cylindrical object), the network was able to anticipate
whether its image would be wider or narrower after a 90° rotation. Blue lines and shaded areas indicate the median and area between the 25th and75th percentiles over the dataset, whereas black lines show trajectories for an object that is shown at the bottom. d, Optimal encoding dimensionality
increasing with training sample size, in the absence (blue) and presence (red) of the large observation noise. e, Comparison of test prediction error,

2  

defined by errork = gt+k − WT ut+k|t  / |gt+k |2 , where gt = st − ωt indicates the observation-noise-free input. PredPCA (solid lines) show a smaller
test prediction error and an earlier error convergence compared with the naive autoregressive model (dashed lines). Blue and red lines denote the error in
the absence and presence of the large observation noise. f, Comparison of test prediction error between PredPCA, naive autoregressive model, TICA, TAE,
DMD and SSMs based on the Kalman filter (KF) and the Bayesian filter (BF), when trained with 800 objects in the absence of noise. d−f are obtained withten different realizations of training and test samples. The green bars in f indicate the minimum test prediction error among these ten different realizations.
The shaded areas and error bars indicate the standard deviation. See Supplementary Methods sections 1 and 2 for further details.



We note that we also trained PredPCA with an image dataset ofrotating 3D human faces and confirmed that PredPCA can accurately predict subsequent images and extract relevant features suchas pan and tilt angles of face images in an unsupervised manner(Supplementary Fig. 2).
As a further application to more natural data, we lastly trained
a neural network with natural scenes captured from a driving car37(Fig. 4a and Supplementary Video 3). Here, we aimed at demonstrating the applicability of our simple, analytically solvable linear
method to real-world video prediction and feature extraction tasks,
rather than comparing the prediction accuracy of PredPCA withthat of state-of-the-art video prediction methods exploiting engineering wisdom. For predictions, we separated the videos into six
groups of data based on the magnitude of change in the images perframe and trained six predictors separately with each group of data;subsequently, post hoc PCA was applied to the synthesized predicted
input (see Supplementary Methods section 1 for further details).
440

PredPCA could predict 0.5-seconds future images of previously
unexperienced natural scenes with a certain accuracy (Fig. 4a).
Moreover, PredPCA could extract brightness, the vertical and lateral asymmetries, and lateral motion in the scenes underlying thedriving car videos (Fig. 4b,c). For the feature extractions, the entire
video was simply supplied to PredPCA without the six sub-groups;
thus, the global convergence was theoretically guaranteed. We
observed tight correspondences between features learned based on
different finite training samples (Fig. 4b,c, insets), implying that
PredPCA could extract features unique to the generative processthat generated sensory data. In particular, the PC1 of dynamical
features that encoded the lateral motion in the scenes (Fig. 4c) was
relevant to predicting the steering of the car. The features extractedusing PredPCA were retained even when using the data grouping(Extended Data Fig. 5a, b). Other major categorical and dynamical features represent different categories of scenes (Extended Data
Fig. 5c) and motions in different positions (Extended Data Fig. 5d),Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCEa

Example 1

2

3

4

Ground truth

Prediction

Latest input

Prediction error

Difference

b

3

0
4
1st half–2
Brightness of scenery

0

3

3

1

0

–3
–0.3

0

2nd half

0
–3
–3

0

1

3

1st half
PC3 of xt

PC2 of xt

0

0.5

–3

1st half

2

0

0
–3

2

0.3

0.6

0

–3
–0.4

Vertical asymmetry0
Lateral asymmetry

0
–2
–2

0

2

1st half
PC1 of ∆xt+3∣t0

Dynamical features3
2nd half

3
2nd half

2nd half

4

PC1 of xt

cCategorical features
4

0.4

0

–1
–0.01

0

0.01

Lateral motion

Fig. 4 | PredPCA of natural scene videos. a, Four examples of predictions. Top row: ground truth or target image (st+15); that is, 0.5-seconds future image ofthe latest input. Second row: predicted image obtained using PredPCA (WTut+15|t). It should be noted that the blurry edges in the predicted images occurred
primarily because PredPCA predicted the mean of future outcome images. Third row: latest input image (st). Fourth row: prediction error between groundtruth and predicted images; the white regions indicate the extent of errors (magnitude of st+15 − WTut+15|t). Bottom row: difference between ground truth andlatest input images (magnitude of st+15 − st). The test prediction error, defined by error15 =

 


st+15 − WT ut+15|t 2 / |st+15 − st |2 , was 0.648. Althoughunexpected events during 0.5 s were unpredictable and some predictions were inaccurate owing to the limited effective dimensionality of the input, theresults indicate that PredPCA provides predictions that interpolate unseen future images and the latest input images, without using any label for training.
b, Extraction of brightness and the vertical and lateral asymmetries in driving car videos as the PC1–PC3 of categorical features (that is, x̄t). Insets depict
tight correspondences between features extracted using PredPCA, learned exclusively on the basis of the first and second half of training samples. For
PC1−PC3, the correlations between them are larger than 0.9999. c, Extraction of lateral motion from driving car videos as the PC1 of dynamical features(that is, Δxt+3|t). The correlation between features learned based exclusively on the first and second half of training samples is 0.9962. In b and c, blue lines
and shaded areas indicate the median and area between the 25th and 75th percentiles. Refer to Supplementary Methods section 1 for further details.

respectively. Remarkably, unlike conventional video prediction
methods6, PredPCA could extract these relevant features in an
unsupervised manner, without the use of labels or target signals
for training.

In summary, using examples of objects rotating in 3D and natural scenes, we demonstrated that PredPCA can filter out observation noise and minimize test prediction error by extracting features
relevant to generalizing predictions. Although the true generative

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

441

Articles

NATurE MACHInE InTEllIGEnCE

process is unknown for these examples, these results indicate thatthe outcomes of PredPCA capture the plausible properties of natural data. These results highlight the prediction generalization andfeature extraction capabilities of PredPCA as well as its wide applicability to real-world data.

Discussion

Our proposed scheme, PredPCA, is thus shown to identify a concise
representation that provides the global minimum of the test prediction error, by first predicting subsequent observations and thenperforming post hoc PCA of the predicted inputs. This is essential
for maximizing the prediction generalization capability, as well as
for ensuring accurate and unbiased estimation of system properties,
comprising hidden states, system parameters and dimensionalities.
Our scheme is formally based on Akaike’s statistics19,27 and is consistent with existing information-theoretical views of biological optimizations, including maximum negentropy38, predictive coding2,3,
predictive information39 and the free energy principle40—providing
a normative, analytically solvable example of a neural network thatmaximizes information quality and generalization capability.
PredPCA offers an interpretable hidden state representation(Methods section ‘Asymptotic linearization theorem') that is preferable for generalizing prediction, without using prior knowledge about external systems. To this end, the global convergenceguarantee or convex optimization of PredPCA (Methods sections‘Derivation of PredPCA' and ‘Test prediction error minimization')
is essential, because the representation could otherwise change
depending on the initial conditions, training history or level ofobservation noise, rendering such representation unreliable, andmay overfit to a particular training dataset. PredPCA further guarantees the asymptotic identification of the true system properties witha global convergence guarantee (Methods sections ‘Asymptotic linearization theorem', ‘System parameter identification' and Table 2)
when sensory data are generated from a class of canonical systems (Methods section ‘System'), provided with sufficient training
samples T and sufficiently high-dimensional observations satisfying Ns ≫ Nx ≫ 1. This is remarkable because PredPCA can identify properties of canonical systems even with a limited number oftraining samples, up to a small range of errors that are inversely proportional to the sample and system sizes, as empirically validated
in Fig. 2 (and Extended Data Fig. 1 and Supplementary Fig. 1; seealso ref. 34). These guarantees are crucial, particularly when feature
extraction failures or misunderstanding of a system lead to catastrophic problems in subsequent applications, such as in automateddriving or medical diagnosis. In general, finite but sufficiently large
T, Ns/Nx, and Nx are required to ensure this asymptotic property—
because only under such conditions are the major principal components of the de-noised input guaranteed to match the hidden statesof the original nonlinear system34.
Unlike PredPCA, conventional (nonlinear) prediction strategies
using autoencoders10,11, TAE13 or SSMs16–18 do not have such guarantees and can fail to reliably provide accurate prediction and feature
extraction depending on various conditions, as shown in Figs. 2 and
3 and Methods section ‘Filtering out observation noise', because
they have many spurious solutions. Having said this, if a learner hasa sufficient amount of prior knowledge about the generative processthat generates sensory data (for example, knowledge about underlying physics), incorporating such knowledge into prediction can provide more interpretable and accurate predictions. Such knowledge
may remove spurious solutions and make all solutions the global
minimum. In other words, for these related methods, the outcomesof PredPCA are potentially of great importance in setting a plausible
initial condition and an appropriate empirical prior, in the absenceof prior knowledge.
It should be noted that if the number of training samples issufficient and the magnitude of observation noise is sufficiently
442

low, the prediction error of PredPCA may be larger than that ofstate-of-the-art prediction methods using deep neural networks4–6—
because the generalization error may be negligible under such a
condition. The improvement in generalization capability obtainedby omitting minor eigenmodes has been reported using deep neural networks41,42. This implies a potential extension of PredPCA toanalytically solve the optimal representations for deep learning in a
weakly nonlinear regime.
Although this work focuses on discrete-time systems, one may
think of these systems as approximations of the physical reality incontinuous time, where generative processes exhibit a hierarchyof timescales. From this perspective, the definitions of signal andobservation noise will change depending on the time bin size ofobservations. Thus, it is crucial for accurate predictions to determine the time bin size to ensure that the timescale of observations
matches the timescale of generative processes.
PredPCA-type learning can be implemented in biological neuronal networks and biologically inspired neuromorphic hardware43,44.
Neurons in these systems must update their synapses to perform predictions under physiological or physical constraints—in particular, it
is difficult for them to access non-local information, such as the synaptic weights of other neurons45. This fact limits biologically plausible
learning to a local rule that updates synapses based on and only on
pre- and post-synaptic neural activities and additional directly accessible signals. Conventional PCA and ICA algorithms are non-local,
but we have previously developed a local learning algorithm that performs both PCA and ICA46,47. This algorithm can make PredPCA alocal learning rule and guarantees its biological plausibility. Hence,
we can speculate that PredPCA-type learning underlies the generalization capability of biological organisms and the self-organization ofinternal models48. Some neuronal substrate, such as neuromodulators49,50, may encode the test prediction error expectation for mediating structural learning or model selection in the brain.
For further discussion, please refer to the Supplementary
Discussion.
In summary, PredPCA has proved to be an analytically solvable
unsupervised dimensionality reduction scheme capable of extracting the most informative components for generalizing prediction.
By effectively filtering out unpredictable noise, PredPCA can reliably identify plausible system properties, with a global convergenceguarantee, and can globally minimize the test prediction error.
Although this paper focuses on the autoregression, PredPCA canminimize the generalization error for a class of regression tasks,
indicating its potential applicability to various real-world applications. As a mathematically proved optimal generalization strategy,
our scheme is potentially useful for understanding biological generalization mechanisms and for creating reliable and explainableartificial general intelligence.

Methods

In what follows, we mathematically express the benefits of PredPCA. Methods
sections ‘System' and ‘Derivation of PredPCA' formally define the system and
PredPCA. Methods sections ‘Filtering out observation noise' and ‘Test prediction
error minimization' prove that PredPCA inherits preferable properties of boththe standard PCA and autoregressive models, and outperforms naive PCA andautoregressive models in terms of robustness to noise and generalization ofprediction. Methods sections ‘Asymptotic linearization theorem' and ‘Systemparameter identification' demonstrate that PredPCA identifies the optimal hidden
state estimator and the true system parameters of a class of canonical systems witha global convergence guarantee, owing to the asymptotic property of linear neuralnetworks with high-dimensional inputs34. Supplementary Methods sections 1 and 2
provide the simulation protocols.
System. We suppose that a system in the external milieu is expressed as xt+1 = ft + zt
and st = gt + ωt, where ft ≡ f(xt, xt-1,…) and gt ≡ g(xt) are nonlinear functions of
xt, while zt and ωt are mutually independent white noises characterized withzero means and covariances Σz and Σω. We assume that the system is in a steadystate. To generate predictable dynamics, Σz is assumed to be smaller than Σx in
magnitude; whereas we typically consider a large Σω. Without loss of generality, we

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCE
can suppose that the steady state of xt follows a distribution with zero mean and theidentity covariance Σx ≡ I. For analysis, we consider a family of functions ft ≡ Bψt
and gt ≡ Aψt spanned by nonlinear basis functions ψ t ≡ ψ (xt ) ∈ RNψ , where Nψdenotes the number of linearly independent bases, B ∈ RNx ×Nψ is a full-row-rank
transition matrix, and A ∈ RNs ×Nψ is a full-column-rank mapping matrix from thebases to the sensory input. Thus, equation (1) becomes
st = Aψ t + ωt

(8)

xt+1 = Bψ t + zt

(9)and equation (2) becomes

As the dimensionality of bases increases, each element of f(xt) and g(xt)asymptotically expresses an arbitrary nonlinear mapping if A and B are suitably
selected (refer to universality). We assume Nx ≤ Nψ ≤ Ns such that the system
dynamics are produced by hidden states that are lower-dimensional than theobservations. Although this paper supposes ψt = ψ(xt), the same analysis can beapplied to a system comprising ψt = ψ(xt, xt-1,…) by redefining (xt, xt-1,…) and (st, st1,…) as new xt and st, respectively. Table 1 presents the glossary of expressions.
Derivation of PredPCA. PredPCA aims to minimize the multistep prediction
error for predicting a 1-to-Kf-step future of the aforementioned system by
optimizing synaptic weight matrices using and only using the current and pastobservations st , st−1 , …, st−Kp +1, where Kf and Kp are imposed by the problem
setup. Hidden states and bases (xt, ψt), system parameters (A, B, Σx, Σψ, Σz, Σω), andthe numbers of hidden state and basis dimensions (Nx, Nψ) are unknown
to a learner.
The error for predicting the k-step future is defined by εt+k|t ≡ st+k − WTVkϕt,()T
where ϕt ≡ sTt , sTt−1 , …, sTt−Kp +1 ∈ RNϕ is a vector of observations,

W ∈ RNu ×Ns is the transpose of the decoding synaptic weight matrix, andVk ∈ RNu ×Nϕ is the kth encoding synaptic weight matrix. Although general
nonlinear bases can be used as ϕt, a simple vector of observations serves the
purpose of this paper. We will show below that the prediction and system
identification using these linear bases are accurate when the dimensionality of
inputs are sufficiently large. Minimizing εt+k|t can be viewed as a generalizationof the standard PCA29 that minimizes the reconstruction error of the current
observation (that is, ϵPCA
≡ st − WT Wst ).
t
Formally, the cost function of PredPCA for multistep predictions is defined by
L≡

Kf 
 

1
ϵt+k|t 2 ,
2
q

(10)

k= 1

∑T
where ⟨•⟩q ≡ T1
t=1 • is the expectation over the empirical distribution q.
Solving the fixed point of the above cost function L with respect to Vk yields the
optimal estimator. From∂L
T
= −W ϵt+k|t ϕt∂Vk⟨minimization of L is achieved by maximizing the second term under the constraintof WWT = I (note that this constraint is automatically satisfied by minimizing L).
Hence, the optimal W is provided as the transpose of the major eigenvectors of. This solution is unique up to the multiplication of an Nu × Nu orthogonalΣPredsmatrix from the left. The global convergence and absence of spurious solutions
are guaranteed even when W is computed by equation (13) because of the global
convergence property of Oja’s subspace rule for PCA30. In short, PredPCA is aconvex optimization and thus it can reliably identify the optimal synaptic weightmatrices W and V1 , …, VKf for predictions, which provides the global minimumof the cost function L.
Filtering out observation noise. Here, we compare the components extractedusing PredPCA and the standard PCA28,29. We show that only PredPCA can remove
observation noise and accurately estimate the observation matrix A as training
sample size T increases.
We introduce the∫ expectation over true distribution p (ϕt , st+1 , …, st+Kf ),denoted by ⟨•⟩ ≡ • p (ϕt , st+1 , …, st+Kf ) dϕt dst+1 · · · dst+Kf . The empirical
distribution approaches this true distribution in the large training sample
size limit: p (ϕt , st+1 , …, st+Kf ) = plim q (ϕt , st+1 , …, st+Kf ). ThroughoutT→∞the manuscript, we suppose 〈st〉 = 0, 〈ψt〉 = 0, and 〈xt〉 = 0 for the sake ofsimplicity. The true
matrix⟨ covariance⟩⟨ ⟩ of some variable ξt is denoted byΣξ ≡ cov [ξt ] ≡ ξt ξTt − ⟨ξt ⟩ ξTt . Here, any estimator or statistic θ underconsideration, calculated based on the empirical distribution, can be decomposed
into its true value θ and its generalization error δθ ≡ θ − θ, where δθ is in the T−1/2
order (see Supplementary Methods section 4 for the conditions and the proof).
Below, we will decompose θ into θ and δθ and then solve θ analytically.
The standard PCA conducts the eigenvalue decomposition of⟨the actual
input⟩
covariance, calculated based on the empirical distribution: Σs ≡ st sTt q. The
convergence to some unknown underlying distribution in the large-sample limit is51. From equation (8),is decomposed asa known property( of PCA)( the covariance)Σs = Σs + O T−1/2 = AΣψ A + Σω + O T−1/2 owing to the independenceof ψt and ωt. As the observation noise covariance Σω is involved in Σs, the major
eigenvectors of Σs that PCA extracts are biased toward the directions of the noise’s
major eigenvectors. This bias is a common issue of autoencoding approaches10,11
that renders the identification of the true system parameters difficult.
In contrast with the standard PCA, PredPCA conducts the eigenvalue⟩
Kf ⟨∑decomposition of the predicted input covariance: ΣPred
st+k|t sTt+k|t .
≡ K1
s
f

where ΣψPred ≡ K1

f

⟩

q

= −W

⟨(

st+k − WT Vk ϕt ϕTt
)

⟩

q

= O,

(11)

under an assumption of WWT = I (which is preserved by equation (13) below), theoptimal Vk is found to be
⟩ ⟨
⟩−1
⟨
ϕt ϕTtVk = W st+k ϕTt
(12)
q

q

Ẇ ∝ −

∂L
=
∂W

Kf ⟨
∑

(
ut+k|t st+k − WT ut+k|t

k=1

T

) ⟩
q

(13)

Equation (13) is similar to Oja’s subspace rule for PCA28 except that st+k|t isused instead of st+k to define ut+k|t. In this sense, PredPCA conducts post hoc
dimensionality reduction (PCA) of the predicted input. The update byequation (13) maintains W as an orthogonal matrix (that is, WWT = I) throughout
the learning.
The above PredPCA solution can also be obtained by eigenvalue
decomposition. When WWT = I, the cost function is transformed as



Kf 
 

st+k − WT Wst+k|t 2 = Kf tr [Σs ] − tr WΣPred WT , where
L = 12s2q⟩Kf ⟨⟨k=1 ⟩∑Σs ≡ st sTt q and ΣPred
st+k|t sTt+k|t are the actual and predicted input
≡ K1f
s
k=1

q

covariances calculated based on the empirical distribution, respectively. Thus, the

q

Kf ⟨
∑

k= 1

⟩
⟨
⟩
ψ t+k ϕTt Σϕ−1 ϕt ψ Tt+k is the predicted hidden

basis covariance, calculated based on the true distribution. Applying theprovides the set of major eigenvectors
eigenvalue
decomposition
to ΣPreds()
Ps ≡ P·1 , …, P·Nψ ∈ RNs ×Nψ that correspond to asymptotically non-zero
eigenvalues of the predicted input covariance. Because of the uniqueness of theeigenvalue decomposition, Ps converges to matrix A as the number of training
samples increases—up to the multiplication of a full-rank matrix Ωψ ∈ RNψ ×Nψ
from the right-hand side. Hence, we refer to Ps as the estimator of A:
A ≡ PsWe define the maximum likelihood estimator of st+k as st+k|t ≡ Qkϕt, where
⟨
⟩ ⟨
⟩ −1
Qk ≡ st+k ϕTt q ϕt ϕTt q is the optimal (maximum likelihood) matrix estimator.

Throughout the Article, a bold variable (for example, st+k|t) indicates the estimatorof the corresponding italic variable (for example, st+k). The kth encoder ut+k|t is thus
defined by ut+k|t ≡ Wst+k|t. The optimal W is determined by the gradient descent
on L:

k=1

Owing to this construction, the identification of system parameters (A, B, Σx,Σψ, Σω, Σz) based on PredPCA is not biased by thenoise. From⟩
( the)⟨ observation
independence between ωt+k and ϕt, st+k|t = A ψ t+k ϕTt Σϕ−1 ϕt + O T−1/2
holds. Thus, we obtain
( 1)
ΣPred
= AΣψPred AT + O T− 2 ,
(14)
s

(

1

= Ps + O T− 2
−1

= AΩψ

(

)

+O T

(15)

−1

2

)

Here, we introduced the inverse of Ωψ (instead of Ωψ itself) for
our convenience. Note that Ps is the set of major eigenvectors of thegeneralization-error-free predicted input covariance ΣsPred ≡ AΣψPred AT. In
short, PredPCA can identify matrix A with asymptotically zero error without
directly observing ψt for large T. Notably, the number of basis dimensions Nψ isalso identifiable by counting the number of asymptotically non-zero eigenvalues, which converges to the true Nψ of canonical systems for a large trainingof ΣPredssample size (see Methods section ‘Test prediction error minimization' for theformal definition of the estimator Nψ using the test prediction error).
In addition, multiplying PTs by the predicted input yields the predicted basis
estimator:
ψ t+k|t ≡ PTs st+k|t(1= Ωψ ψ t+k ϕTt Σϕ−1 ϕt + O T− 2

⟨

⟩

)

(16)

The last equality holds from the orthogonality of eigenvectors, that is,
PTs A = PTs Ps Ωψ = Ωψ , and the independence between ωt+k and ϕt. Indeed, ut+k|tNature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

443

Articles

NATurE MACHInE InTEllIGEnCE

with optimized synaptic weight matrices is equivalent to ψt+k|t when Nu = Nψ. In
short, PredPCA can provide the maximum likelihood estimator of the hidden bases
without directly observing ψt+k —up to the multiplication of the full-rank ambiguity
factor Ωψ from the left-hand side. This ambiguity factor is safely absorbed into thedefinition of ψt, without changing the system dynamics, by applying the following
transformations: Ωψψt → ψt, Ps = AΩψ−1 → A, and BΩψ−1 → B. Therefore, theestimated hidden dynamics are formally homologous to the original dynamics.
In terms of conceptual innovations of PredPCA, our analyses reveal that thisscheme can identify the true hidden states, parameters, and dimensionalities ofa class of canonical systems (see below). In particular, the multi-time-step bases
function ϕt is an essential difference between PredPCA and related methods suchas TICA14, TAE13 and DMD15,52. Empirical observations highlight the importanceof filtering out observation noise to reliably perform system identification (Fig. 2,
and Extended Data Figs. 1 and 2). Indeed, features extracted from TICA or DMD
are expressed as complex numbers, which do not match the true hidden states.
Although TAE can identify matrix A and the extracted features are denoted in real
numbers, it still fails to identify true hidden states and other parameters because it
fails to filter out large observation noise (Fig. 2b and Extended Data Fig. 2a).
Furthermore, we presented an algorithm to update synaptic weights(equation (5)), which makes it easier to design a computational architecture for
PredPCA. As discussed above, it is fairly straightforward to implement PredPCA
in neuromorphic hardware through a previously developed local learningalgorithm46,47,53—wherein a previous work has implemented the local algorithm
using resistive random-access memories44. It should be emphasized that PredPCA
is suitable for neuromorphic hardware relative to TCIA, TAE and DMD because
the computations for inverse matrices, complex numbers, and eigenvalue
decomposition of non-symmetric matrices are intractable in neural networks. For
further discussion, please refer to Supplementary Discussion.
Test prediction error minimization. A learner needs to predict the futureconsequences of unseen input data based on learning with a limited number oftraining samples. Here, we analytically solve the expectation of the PredPCA’s test
prediction error as a function of the training samples (T), encoding dimensions (Ns),and number of past observations used for prediction (Nϕ = KpNs). Its minimization
enables a learner to maximize the generalization ability by optimizing free parametersin the network without knowing the true distribution that generates test samples.
PredPCA’s test prediction error is defined as the squared error over the true
distribution p. Meanwhile, the learning is based on the empirical distribution q.
Thus, the test error is given as a functional of q:
Ltest [q] ≡

K

f 
 

1
ϵt+k|t [q]2
2

(17)

k=1

Here, the prediction error (which is also a functional of q) is given
as ϵt+k|t [q] ≡ st+k − Ps PTs st+k|t using the major eigenvectors of
the predicted input covariance Ps ≡ (P·1 , …, P·Nu ) ∈ RNs ×Nu and the
⟨
⟩ ⟨
⟩ −1
maximum likelihood estimator st+k|t = st+k ϕTt q ϕt ϕTt q ϕt computed

based on the empirical distribution q. The generalization error of major
eigenvectors Ps is negligible up to the leading order (see Supplementary Methods
section 5 for details). The q-dependent factor in st+k|t is computed as⟨⟩ ⟨⟩ −1 ( ⟨⟩⟨⟩ )(⟨⟩ ) −1st+k ϕTt q ϕt ϕTt q = st+k ϕTt + δ st+k ϕTt q
Σϕ + δ ϕt ϕTt q⟨() ⟩= Qk + δ st+k − Qk ϕt ϕTt q Σϕ−1 up to the leading⟨⟩
order, using the optimal mapping Qk ≡ st+k ϕTt Σϕ−1 (notethat δ〈∙〉q ≡ 〈∙〉q − 〈∙〉). Thus, the prediction error becomes) ⟩⟨(ϵt+k|t [q] = st+k − Ps PTs Qk ϕt − Ps PTs δ st+k − Qk ϕt ϕTt q Σϕ−1 ϕt ,where Ps ≡ (P·1 , …, P·Nu ) ∈ RNs ×Nu denotes the major eigenvectors of thegeneralization-error-free predicted input covariance ΣsPred. Then, we define theexpectation of Ltest[q] over different empirical distributions q, given by
L ≡ E{q} [Ltest [q]](18)Here, E{q}[⋅] means the expectation over different empirical distributions. The
expectation over different q is equivalent to the[expectation
over )p for⟩a]linear term⟨(that involves a single δ〈∙〉q factor. Hence, E{q} δ st+k − Qk ϕt ϕTt q = O.
In contrast, a term that comprises the square of δ〈∙〉q yields the positive variancethrough[the interaction of the two factors, which is computed ]
⟨(
) ⟩
⟨(
) ⟩T
as E{q} δ st+k − Qk ϕt ϕTt q Σϕ−1 δ st+k − Qk ϕt ϕTt q)
N (= Tϕ Σs − Qk Σϕ QTk . Therefore, we find


Kf 
L
=
tr [Σs ] − tr PTs ΣsPred Ps

2



test error expectation
training error

 
 3
Kf N ϕ  T 
tr Ps Σs − ΣsPred Ps + O T− 2
+
2T



generalization error

444

(19)

See Supplementary Methods section 5 for the detailed derivation. Thisis viewed as a variant of AIC19 and NIC20. For practical use, covariances andeigenvectors in equation (19) are replaced with their estimators: Σs → Σs,, and Ps → Ps, whereΣPred→ ΣPredss[ ( L does not change) ] by these replacements inthe leading order. Because tr PTs Σs − ΣsPred Ps > 0, the generalization error

monotonically increases with the dimensionality of the encoders Nu. Meanwhile,
the reduction of the training prediction error becomes small as Nu increases, andit reaches zero for Nu > Nψ due to zero eigenvalues of ΣsPred. Hence, the optimal Nu
that minimizes L is less than Ns.
The optimal encoding dimensionality that minimizes L is comparable to theeffective dimensionality of true hidden basis dynamics of canonical systems forlarge T. Thus, Nψ ≡ argminNu L provides the estimator of the true hidden basis
dimensionality. In particular,] when T is larger than a large[ Nψ = Nψ holds

finite constant Tcψ ≡ Nϕ tr Σs − ΣsPred / (Λs )Nψ Nψ , where (Λs )Nψ Nψ is the Nψth(that is, the smallest non-zero) eigenvalue of ΣsPred. In contrast, equation (19)(AR) modelwith Nu = Ns provides the test prediction error of(an autoregressive
) [
]
K
N
that does not consider hidden states: LAR = 2f 1 + Tϕ tr Σs − ΣsPred .

Because[ ( some components) ] of Σ[ω are generally] perpendicular to Ps,
tr PTs Σs − ΣsPred Ps < tr Σs − ΣsPred for Nu < Ns. This means that the test
prediction error of PredPCA with optimal Nu is smaller than that of autoregressivemodels. Hereafter, we suppose Nu = Nψ = Nψ.

Asymptotic linearization theorem. The asymptotic linearization theorem34
was originally introduced to guarantee accurate extraction of independently
and identically distributed hidden sources from its high-dimensional nonlinear
transformations. In this Article, we use this theorem to prove that the true hidden
state xt ∈ RNx is accurately estimated from its unknown nonlinear transformation
ψ (xt ) ∈ RNψ with asymptotically zero element-wise error as Nx and Nψ/Nx (andT) diverge. In this section, we suppose that ψ(xt) is expressed in a specific but
generic form of two-layered structure, ψ(xt) = Cρ(Rxt + r). Here, the elements ofR ∈ RNψ ×Nx and r ∈ RNψ are fixed Gaussian random variables independently
drawn from N [0, 1/Nx ]; C ∈ RNψ ×Nψ is a matrix whose elements are, on
1/2
average, on the order of N−; and ρ (·) : R �→ R is an odd nonlinear function,
ψwhere the correlation between ρ(ξ) and a unit Gaussian variable ξ sampledfrom p(ξ) ≡ N [0, 1] is not close to zero. When Nψ is large, each elementof ψ(xt) can represent an arbitrary nonlinear mapping of xt by adjusting C(refer to universality)54–57.
The assumption behind the theorem is as follows: (1) the elements of hidden
states xt are not strongly dependent on each other (where zero mean and identity
covariance matrix are supposed without loss of generality), in the sense that the
average of higher-order correlations of the elements of xt asymptotically vanishes
for large Nx with less than the order of 1; and that (2) the matrix components of
C that are parallel to R are not too small compared to the other components (thatis, the mapping is not very close to a singular mapping)—namely, the ratio of theminimum eigenvalue of RTCTCR to the maximum eigenvalue of CCT is assumed to
be much greater than 1. We note that RT R = O (Nψ /Nx ) is much greater than 1,
so the condition (2) is easily satisfied when singular values of C are of order 1. The
asymptotic linearization theorem states that under these two conditions, covarianceΣψ has a clear spectrum gap that separates major and minor components, where themajor and minor components correspond to linear and nonlinear transformationsof the true hidden states, respectively.
Let P ∈ RNψ ×Nx be the set of the first to Nxth major eigenvectorsof Σψ, and let Λ ∈ RNx ×Nx be the diagonal matrix of the corresponding
eigenvalues. The asymptotic linearization theorem proved that applying
PCA to ψ(xt) provides accurate estimation of xt up to the multiplication of
a fixed√orthogonal matrix Ω; that is, Λ−1/2 PT ψ (xt ) = Ω xt + O (σ x ). Here,
(
)
(
)
2
2
2
σx =
ρ2 /ρ′ − 1 (1 + λ) Nx /Nψ + ρ′′′ / 2ρ′ Nx is the standard∫
deviation of the linearization error, where ρ2 ≡ ρ2 (ξ) p (ξ) dξ,∫ dρ(ξ)∫ d3 ρ(ξ)′′′ ≡are statistics of the nonlinearρ′ ≡dξ p (ξ ) dξ , and ρdξ 3 p (ξ ) dξfunction ρ over unit Gaussian variable ξ, and λ is an order-one constantdetermined by the characteristics of C. The linearization error monotonicallydecreases as the system size increases (that is, when Nψ/Nx and Nx diverge)—
asymptotically achieving the zero-element-wise-error hidden state estimation in
the large system size limit. Please refer to ref. 34 for further details.
This theorem can be applied to the estimator of ψ(xt). Let Pψ ∈ RNψ ×Nx be themajor eigenvectors of Σψ (see equation (23) below for its definition and analyticalsolution), and Λψ ∈ RNx ×Nx be the corresponding eigenvalues. The hidden state
estimator is given as−1

xt+k|t ≡ Λψ 2 PTψ ψ t+k|t
−1

(

1

= Λψ 2 PTψ Ωψ ψ t+k ϕt Σϕ−1 ϕt + O T− 2

⟨

⟩

)

(20)

From the asymptotic linearization theorem,

1/2 T
Λ−
Pψ Ωψ ψ t+k = Ωx xt+k + O (σ x ) holds, where Ωx ∈ RNx ×Nx is a fixedψ

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCE
orthogonal matrix. Here, we treated Ωψψt+k as a nonlinear function of xt+k andapplied the theorem. Thus, equation (20) is solved analytically as( 1)⟩⟨
xt+k|t = Ωx xt+k ϕt Σϕ−1 ϕt + O T− 2 + O (σ x )(21)This⟩result shows that the maximum likelihood estimator of xt+k based on ϕt,
xt+k ϕt Σϕ−1 ϕt , is available (up to the ambiguity factor Ωx, and the order T−1/2
and σx small error terms), despite the fact that PredPCA is unsupervised learningthat does not use any explicit data of xt+k for training. Similar to Ωψ, the ambiguityof Ωx can be absorbed into the definition of xt, without changing any system
dynamics, by applying the following transformations: Ωxxt → xt, ΩxB → B, Ωxzt →zt, and RΩx−1 → R. Notably, the number of state dimensions Nx is also identifiable
by defining the estimator Nx as the largest spectrum gap of Σψ, which is guaranteed
to converge to true Nx when σx is smaller than a small positive constant σ cx and T islarger than a large finite constant Tcx .
It is well known that conventional nonlinear blind source separation
approaches using nonlinear neural networks (for example, nonlinear ICA) do notguarantee the identification of true hidden sources under the general nonlinear
blind source separation setup58,59. In contrast, it is remarkable that the asymptotic
linearization theorem mathematically guarantees the achievability of the nonlinear
blind source separation when Nψ ≫ Nx ≫ 1 (ref. 34).

⟨

System parameter identification. We demonstrated above that PredPCA canidentify the true observation matrix A. Here, we show that it can also identify othersystem parameters B, Σψ, Σx, Σω and Σz asymptotically—if the assumptions of theasymptotic linearization theorem are met and the number of training samples islarge.
These parameter identifications are based on the linearized transition mappingfrom ψt to ψt+1, denoted by Ψ; thus, we first compute the estimator ⟨of Ψ. We⟩
decompose ψt+1 as ψ t+1 = Ψ ψ t + Δψ t+1|t + O (zt ), where Ψ = ψ t+1 ψ Tt Σψ−1is the optimal basis transition matrix, Δψt+1|t is the linearization error that isperpendicular to both ψt and zt, and O (zt ) is a term related to small noise zt. ThisΨ can be viewed as a finite basis size version of the Koopman operator52,60. The
Tbasis estimator based on the current input( is defined) as ψ t|t ≡ Ps st andcomputed as ψ t|t = Ωψ ψ t + PTs ωt + O T−1/2 . Using this, we have
⟩
⟨(
⟨
)(
)T ⟩
(
)
+ O T−1/2
ψ t+k|t+k ψ Tt|t = Ωψ ψ t+k + PTs ωt+k Ωψ ψ t + PTs ωt
q
⟨
⟩
(
)
noise
is white= Ωψ ψ t+k ψ Tt ΩψT + O T−1/2 as the observation⟨⟩and independent of ψt andψ
.
Inparticular,
ψ t+1 ψ Tt = Ψ Σψ and
t+k
⟨
⟩
⟨
⟩
ψ t+2 ψ Tt = Ψ 2 Σψ + Δψ t+2|t+1 ψ Tt hold. Thus, we obtain the following
estimator of the basis transition matrix:
⟨
⟩ ⟨
⟩− 1
ψ t+1|t+1 ψ Tt|t
Ψ ≡ ψ t+2|t+2 ψ Tt|tq−1

= Ωψ Ψ Ωψ

−1= Ωψ Ψ Ωψ

+ Ωψ

⟨

(

q

Δψ t+2|t+1 ψ Tt

+O T

−

)
1

⟩

(

1

Σψ−1 Ψ −1 Ωψ−1 + O T− 2

)

(22)

+ O (σ ψ )

2

⟨
⟩
⟨
⟩
zt+2 and ϕ(t, xt+2 ϕTt ) = B⟨ ψ t+1 ϕTt⟩ holds. Thus,(equation) (21) for k = 2 becomes−1
xt+2|t = Ωx BΩψ Ωψ ψ t+1 ϕt Σϕ−1 ϕt + O T−1/2 + O (σ x ). Hence, using
equation (16), we find the following estimator of the transition matrix:
⟨
⟩ ⟨
⟩−1
ψ t+1|t ψ Tt+1|t
B ≡ xt+2|t ψ Tt+1|t
q
q
(25)
( 1)
−
−1
= Ωx BΩψ + O T 2 + O (σ x )The hidden state covariance estimator is given by Σx ≡ Σx ≡ I as we definedΣx so. Thus, as equation (9) yields Σx = BΣψBT + Σz, the process noise covarianceestimator is given byΣz ≡ Σx − BΣ ψ BT(1

= Ωx Σz ΩxT + O T− 2)

+ O (σ x )(26)

In summary, PredPCA could identify the true system parameters A, B, Ψ,Σψ, Σx, Σω and Σz with a global convergence guarantee as the system and training
sample sizes increase, using noisy observations only—up to the full-rank linear
transformations (Ωψ, Ωx) that do not change the system dynamics. When zt and ωt
are white Gaussian noises, these parameters are sufficient to identify the canonicalsystem. The aforementioned analyses hold true even when zt and ωt are white
non-Gaussian noises, although in this case, unsupervised identification of thethird- or higher-order moments of zt and ωt has not yet been established. The
zero-element-wise-error identification of these parameters will be asymptoticallyachieved when Nψ/Nx, Nx and T diverge. This global convergence guarantee isan advantage of PredPCA compared with conventional system identification
approaches18,61. If ψ(xt) is a linear function of xt, the true system becomes a linear
system and thus provides σx = σψ = 0; hence, PredPCA is guaranteed to identify
the true system parameters with zero error as the increasing training samples,when Nx ≤ Ns. Table 2 summarizes the definitions and analytical solutions of theseestimators. Every estimator can be computed by following the definition, where its
analytical solution and accuracy have been proved theoretically.
The identification of system properties using PredPCA was empiricallydemonstrated with the example of handwritten digit images (Fig. 2, Extended Data
Fig. 1 and Supplementary Fig. 1). Although it is difficult to prove what the true
generative process is for rotating 3D objects (Fig. 3 and Supplementary Fig. 2)
or natural scenes (Fig. 4 and Extended Data Fig. 5), empirical observations
suggest that PredPCA can extract features relevant to generalized predictions.
At least, PredPCA was able to identify the angles of rotating objects (Fig. 3c andSupplementary Fig. 2) and lateral motion in natural scenes (Fig. 4c and Extended
Data Fig. 5b), indicating the identification of a part of their generative processes.
These observations imply that the outcomes of PredPCA capture the plausible
properties of natural data.
Reporting Summary. Further information on research design is available in the
Nature Research Reporting Summary linked to this article.

Data availability

This estimator converges to the optimal Ψ up to the ambiguity of Ωψ for large
T and Nψ. The variance of the linearization error term O (σ ψ ) is of the same orderas the variance of nonlinearly transformed components of xt that are involved in ψt;34thus, using the asymptotic linearization theorem
the variance of√( , we compute)ρ2 − ρ′2 /Nψ as the orderthe nonlinear components and obtain σ ψ =Image data used in this work are available in the MNIST dataset33 (http://yann.
lecun.com/exdb/mnist/index.html, for Fig. 2), the ALOI dataset36 (http://aloi.
science.uva.nl, for Fig. 3), and the BDD100K dataset37 (https://bdd-data.berkeley.
edu, for Fig. 4). Figures 2–4 are generated by applying our scripts (see below) tothese image data.

(see Supplementary Methods section 6 for further details).
Next, we compute the covariance matrices of hidden bases
and observation⟩noise. By multiplying the inverse of Ψ with⟨()
ψ t+1|t+1 ψ Tt|t = Ω ψ Ψ Σψ ΩψT + O T−1/2 , we find the hidden basis

Code availability

MATLAB scripts used in this work are available at https://github.com/
takuyaisomura/predpca or https://doi.org/10.5281/zenodo.4362249. The scripts arecovered under the GNU General Public License v3.0.

qcovariance estimator (symmetrized version) as
)
(
⟩
⟨
⟩
⟨
Σ ψ ≡ 12 Ψ−1 ψ t+1|t+1 ψ Tt|t + ψ t|t ψ Tt+1|t+1 Ψ−Tq= Ωψ Σψ ΩψT + O

(

−1

T

2

)

q

Received: 5 March 2020; Accepted: 27 January 2021;
Published online: 12 April 2021
(23)

+ O (σ ψ )See Supplementary Methods section 6 for the order of the linearization error
term. Using this Σψ, the observation noise covariance estimator is given asΣ ω ≡ Σs − AΣ ψ AT(1

= Σs − AΣψ AT + O T− 2(= Σω + O T

−1
2

)

)

+ O (σ ψ )(24)+ O (σ ψ )Finally, we estimate the state transition matrix and covariance matrices ofhidden states and process noise. From equation (9) and the independence between

References

1. Rao, R. P. & Ballard, D. H. Predictive coding in the visual cortex: a functional
interpretation of some extra-classical receptive-field effects. Nat. Neurosci. 2,
79–87 (1999).
2. Rao, R. P. & Sejnowski, T. J. Predictive sequence learning in recurrent
neocortical circuits. Adv. Neural Info. Proc. Syst. 12, 164–170 (2000).
3. Friston, K. A theory of cortical responses. Phil. Trans. R. Soc. Lond. B 360,
815–836 (2005).
4. Srivastava, N., Mansimov, E. & Salakhudinov, R. Unsupervised learningof video representations using LSTMs. In Int. Conf. Machine Learning
843−852 (ML Research Press, 2015).
5. Mathieu, M., Couprie, C. & LeCun, Y. Deep multi-scale video prediction
beyond mean square error. Preprint at https://arxiv.org/abs/1511.05440(2015).

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

445

Articles

NATurE MACHInE InTEllIGEnCE

6. Lotter, W., Kreiman, G. & Cox, D. Deep predictive coding networks for video
prediction and unsupervised learning. Preprint at https://arxiv.org/abs/
1605.08104 (2016).
7. Hurvich, C. M. & Tsai, C. L. Regression and time series model selection insmall samples. Biometrika 76, 297–307 (1989).
8. Hurvich, C. M. & Tsai, C. L. A corrected Akaike information criterion forvector autoregressive model selection. J. Time Series Anal. 14, 271–279 (1993).
9. Cunningham, J. P. & Ghahramani, Z. Linear dimensionality reduction: survey,
insights, and generalizations. J. Mach. Learn. Res. 16, 2859–2900 (2015).
10. Hinton, G. E. & Salakhutdinov, R. R. Reducing the dimensionality of datawith neural networks. Science 313, 504–507 (2006).
11. Kingma, D. P. & Welling, M. Auto-encoding variational Bayes. Preprint at
https://arxiv.org/abs/1312.6114 (2013).
12. Hochreiter, S. & Schmidhuber, J. Long short-term memory. Neural Comput.
9, 1735–1780 (1997).
13. Wehmeyer, C. & Noé, F. Time-lagged autoencoders: deep learning of slow
collective variables for molecular kinetics. J. Chem. Phys. 148, 241703 (2018).
14. Pérez-Hernández, G., Paul, F., Giorgino, T., De Fabritiis, G. & Noé, F.
Identification of slow molecular order parameters for Markov model
construction. J. Chem. Phys. 139, 015102 (2013).
15. Klus, S. et al. Data-driven model reduction and transfer operator
approximation. J. Nonlinear Sci. 28, 985–1010 (2018).
16. Kalman, R. E. A new approach to linear filtering and prediction problems.
J. Basic Eng. 82, 35–45 (1960).
17. Julier, S. J. & Uhlmann, J. K. New extension of the Kalman filter to nonlinear
systems. In Signal Processing, Sensor Fusion, And Target Recognition VI Vol.
3068, 182−193 (International Society for Optics and Photonics, 1997).
18. Friston, K. J., Trujillo-Barreto, N. & Daunizeau, J. DEM: A variational
treatment of dynamic systems. NeuroImage 41, 849–885 (2008).
19. Akaike, H. A new look at the statistical model identification. IEEE Trans.
Automat. Contr. 19, 716–723 (1974).
20. Murata, N., Yoshizawa, S. & Amari, S. I. Network information criterion—
determining the number of hidden units for an artificial neural networkmodel. IEEE Trans. Neural Netw. 5, 865–872 (1994).
21. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 6, 461–464 (1978).
22. Vapnik, V. Principles of risk minimization for learning theory. Adv. Neural
Info. Proc. Syst. 4, 831–838 (1992).
23. Arlot, S. & Celisse, A. A survey of cross-validation procedures for model
selection. Stat. Surv. 4, 40–79 (2010).
24. Comon, P. & Jutten, C. (eds) Handbook of Blind Source Separation:
Independent Component Analysis And Applications (Academic Press, 2010).
25. Ljung, L. System Identification: Theory for the User 2nd edn(Prentice-Hall, 1999).
26. Schoukens, J. & Ljung, L. Nonlinear system identification: a user-oriented
roadmap. Preprint at https://arxiv.org/abs/1902.00683 (2019).
27. Akaike, H. Prediction and entropy. In Selected Papers of Hirotugu Akaike387−410 (Springer, 1985).
28. Oja, E. Neural networks, principal components, and subspaces. Int. J. NeuralSyst. 1, 61–68 (1989).
29. Xu, L. Least mean square error reconstruction principle for self-organizing
neural-nets. Neural Netw. 6, 627–648 (1993).
30. Chen, T., Hua, Y. & Yan, W. Y. Global convergence of Oja’s subspacealgorithm for principal component extraction. IEEE Trans. Neural Netw. 9,
58–67 (1998).
31. Bell, A. J. & Sejnowski, T. J. An information-maximization approach to blindseparation and blind deconvolution. Neural Comput. 7, 1129–1159 (1995).
32. Amari, S. I., Cichocki, A. & Yang, H. H. A new learning algorithm for blind
signal separation. Adv. Neural Info. Proc. Syst. 8, 757–763 (1996).
33. LeCun, Y., Bottou, L., Bengio, Y. & Haffner, P. Gradient-based learning
applied to document recognition. Proc. IEEE 86, 2278–2324 (1998).
34. Isomura, T. & Toyoizumi, T. On the achievability of blind source separation
for high-dimensional nonlinear source mixtures. Preprint at https://arxiv.org/
abs/1808.00668 (2018).
35. Dimigen, O. Optimizing the ICA-based removal of ocular EEG artifacts fromfree viewing experiments. Neuroimage 207, 116117 (2020).
36. Geusebroek, J. M., Burghouts, G. J. & Smeulders, A. W. The Amsterdam
library of object images. Int. J. Comput. Vis. 61, 103–112 (2005).
37. Yu, F. et al. BDD100K: a diverse driving video database with scalable
annotation tooling. Preprint at https://arxiv.org/abs/1805.04687 (2018).
38. Schrödinger, E. What Is Life? The Physical Aspect of the Living Cell and Mind
(Cambridge Univ. Press, 1944).
39. Palmer, S. E., Marre, O., Berry, M. J. & Bialek, W. Predictive information in asensory population. Proc. Natl Acad. Sci. USA 112, 6908–6913 (2015).
40. Friston, K., Kilner, J. & Harrison, L. A free energy principle for the brain. J.
Physiol. Paris 100, 70–87 (2006).
41. Oymak, S., Fabian, Z., Li, M. & Soltanolkotabi, M. Generalization guarantees
for neural networks via harnessing the low-rank structure of the Jacobian.
Preprint at https://arxiv.org/abs/1906.05392 (2019).

446

42. Suzuki, T. et al. Spectral-pruning: compressing deep neural network via
spectral analysis. Preprint at https://arxiv.org/abs/1808.08558 (2018).
43. Neftci, E. Data and power efficient intelligence with neuromorphic learning
machines. iScience 5, 52–68 (2018).
44. Fouda, M., Neftci, E., Eltawil, A. M. & Kurdahi, F. Independent component
analysis using RRAMs. IEEE Trans. Nanotech. 18, 611–615 (2018).
45. Lee, T. W., Girolami, M., Bell, A. J. & Sejnowski, T. J. A unifying
information-theoretic framework for independent component analysis.
Comput. Math. Appl. 39, 1–21 (2000).
46. Isomura, T. & Toyoizumi, T. A local learning rule for independent
component analysis. Sci. Rep. 6, 28073 (2016).
47. Isomura, T. & Toyoizumi, T. Error-gated Hebbian rule: a local learning rule
for principal and independent component analysis. Sci. Rep. 8, 1835 (2018).
48. Dayan, P., Hinton, G. E., Neal, R. M. & Zemel, R. S. The Helmholtz machine.
Neural Comput. 7, 889–904 (1995).
49. Frémaux, N. & Gerstner, W. Neuromodulated spike-timing-dependent
plasticity, and theory of three-factor learning rules. Front. Neural Circuits 9,
85 (2016).
50. Kuśmierz, Ł., Isomura, T. & Toyoizumi, T. Learning with three factors:
modulating Hebbian plasticity with errors. Curr. Opin. Neurobiol. 46,
170–177 (2017).
51. Zhu, B., Jiao, J. & Tse, D. Deconstructing generative adversarial networks.
IEEE Trans. Inf. Theory 66, 7155–7179 (2020).
52. Lusch, B., Kutz, J. N. & Brunton, S. L. Deep learning for universal linear
embeddings of nonlinear dynamics. Nat. Commun. 9, 4950 (2018).
53. Isomura, T. & Toyoizumi, T. Multi-context blind source separation byerror-gated Hebbian rule. Sci. Rep. 9, 7127 (2019).
54. Hornik, K., Stinchcombe, M. & White, H. Multilayer feedforward networks
are universal approximators. Neural Netw. 2, 359–366 (1989).
55. Barron, A. R. Universal approximation bounds for superpositions of asigmoidal function. IEEE Trans. Info. Theory 39, 930–945 (1993).
56. Rahimi, A. & Recht, B. Uniform approximation of functions with randombases. In Proc. 46th Ann. Allerton Conf. on Communication, Control, andComputing 555−561 (2008).
57. Rahimi, A. & Recht, B. Weighted sums of random kitchen sinks: replacing
minimization with randomization in learning. Adv. Neural Info. Process. Syst.
21, 1313–1320 (2008).
58. Hyvärinen, A. & Pajunen, P. Nonlinear independent component analysis:existence and uniqueness results. Neural Netw. 12, 429–439 (1999).
59. Jutten, C. & Karhunen, J. Advances in blind source separation (BSS) andindependent component analysis (ICA) for nonlinear mixtures. Int. J. NeuralSyst. 14, 267–292 (2004).
60. Koopman, B. O. Hamiltonian systems and transformation in Hilbert space.
Proc. Natl Acad. Sci. USA 17, 315–318 (1931).
61. Ljung, L. Asymptotic behavior of the extended Kalman filter as a parameter
estimator for linear systems. IEEE Trans. Automat. Contr. 24, 36–50 (1979).

Acknowledgements

We are grateful to S.-I. Amari for discussions. This work was supported by RIKEN
Center for Brain Science (T.I. and T.T.), Brain/MINDS from AMED under grant number
JP20dm020700 (T.T.), and JSPS KAKENHI under grant number JP18H05432 (T.T.). The
funders had no role in the study design, data collection and analysis, decision to publish,
or preparation of the manuscript.

Author contributions

T.I. conceived and designed PredPCA, performed the mathematical analyses andsimulations, and wrote the manuscript. T.T. supervised T.I. from the early state of thiswork, confirmed the rigour of the mathematical analyses and wrote the manuscript.

Competing interests

The authors declare no competing interests.

Additional information

Extended data is available for this paper at https://doi.org/10.1038/s42256-021-00306-1.
Supplementary information The online version contains supplementary material
available at https://doi.org/10.1038/s42256-021-00306-1.
Correspondence and requests for materials should be addressed to T.I. or T.T.
Peer review information Nature Machine Intelligence thanks the anonymous reviewers
for their contribution to the peer review of this work.
Reprints and permissions information is available at www.nature.com/reprints.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
© The Author(s), under exclusive licence to Springer Nature Limited 2021,
corrected publication 2021

Nature Machine Intelligence | VOL 3 | May 2021 | 434–446 | www.nature.com/natmachintell

NATurE MACHInE InTEllIGEnCE

Extended Data Fig. 1 | See next page for caption.

Nature Machine Intelligence | www.nature.com/natmachintell

ArticlesArticles

NATurE MACHInE InTEllIGEnCE

Extended Data Fig. 1 | Supplementary results of PredPCA with handwritten digit images. a, Transition mapping estimated using PredPCA B ∈ R10×10
accurately matches the true transition mapping B ∈ R10×10 that generates the ascending order sequence. Elements of xt+1|t are permuted and sign-flippedfor visualization purpose. b, This is also the case for the nonlinear dynamics. The estimated mapping from xt|t-1 ⊗ xt-1|t-2 to xt+1|t, B̃ ∈ R10×100, wasobtained using the outcomes of PredPCA, which accurately matches the true mapping of the Fibonacci sequence B̃ ∈ R10×100. Here, ⊗ indicates the
Kronecker product. These results indicate that PredPCA offers the identification of the transition rules underlying the linear and nonlinear dynamics,
without observing the true hidden states xt. c, Prediction error in the absence of random replacement and/or monochrome inversion of digit images, as acounterpart of Fig. 2d. PredPCA’s outcomes are retained with or without those distortions, and relevant encoders comprise up to 10 dimensions owing to
the construction of the input data, highlighting the robustness of PredPCA to various types of large noise. In particular, in the presence of monochrome
inversion, irrespective of random replacement of digits, Nu = 10 provides the global minimum of both equations (6) and (7). Conversely, in the absence ofmonochrome inversion, Nu = 9 provides their global minimum as in this case, the 10-dimensional hidden state representation becomes redundant. This isbecause without monochrome inversion, true hidden states take only 10 different positions in the 10-dimensional coordinate, which can be fully expressedby the 9-dimensional coordinate. Remarkably, PredPCA could detect their difference. Note that monochrome inversion corresponds to the first principal
component (PC1) of PredPCA. This is because whether the next image is a ‘black digit on white background’ or ‘white digit on black background’ is themost predictable feature as the monochrome inversion rarely occurs. Thus, a relatively large prediction error in the absence of monochrome inversion isdue to the lack of the PC1. d, PredPCA increases its performance as the number of past observations used for prediction (Kp) increases until reaching its
finite optimum. Left panel: error in categorizing digits, which converges to near zero as Kp increases (refer to Fig. 2b). Middle panel: parameter estimation
error (refer to Fig. 2c). Right panel: test prediction error (refer to Fig. 2d). The blue line is the optimal test prediction error computed via supervised
learning. The red line indicates the theoretical value computed using equation (7), wherein Kp = 10 (green line) gives its minimum, which matches
empirical observations (black circles). These observations imply that predicting single-time-step future outcomes (st+1) using multi-time-step pastobservations (ϕt) is key to reducing those errors. Note that an extension of PredPCA for multi-time-step prediction while retaining its accuracy is provided
in Methods section ‘Derivation of PredPCA'. c and d are obtained with 20 different realizations of digit sequences.

Nature Machine Intelligence | www.nature.com/natmachintell

NATurE MACHInE InTEllIGEnCE

ArticlesExtended Data Fig. 2 | Comparison with related methods. The errors in estimating system parameters (left and middle panels, as a counterpart of Fig. 2c)and in predicting one-step future inputs in test ascending sequence (right panels, refer to Fig. 2d) are shown. a, Performance of linear TAE. Although
it estimates matrix A with high accuracy, it fails to estimate other parameters, because linear TAE (same as PredPCA with ϕt = st) does not effectively
filter out observation noise. Moreover, linear TAE yields a larger test prediction error even relative to PredPCA with ϕt = st owing to the difference in theircost functions. This is because PredPCA (even with ϕt = st) extracts components most important to predicting high variant signals preferentially, andthereby provides the global minimum of the squared error in predicting the non-normalized target signal (under the constraint of ϕt = st), while linear
TAE minimizes a normalized target signal (see Methods section ‘Filtering out observation noise' for more details). For reference, the blue and red lines
in the right panel represent the optimal test prediction error computed via supervised learning and that of PredPCA with ϕt = st, respectively. The results
are obtained with 20 different realizations of digit sequences. b, Performance of SSM based on Kalman filter. SSM also tends to fail system identification
depending on initial conditions and training history, which leads to a relatively larger prediction error. In the left panel, lines and shaded areas indicate themedian and the 25th to 75th percentile area, respectively. The results are obtained with 100 different realizations of digit sequences.

Nature Machine Intelligence | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCE

Extended Data Fig. 3 | See next page for caption.

Nature Machine Intelligence | www.nature.com/natmachintell

NATurE MACHInE InTEllIGEnCE

ArticlesExtended Data Fig. 3 | Accuracy of long-term predictions. PredPCA and SSM can both yield generative models to predict an arbitrary future. However,
SSM can fail to identify system parameters depending on initial conditions and training history, leading to the failure of long-term predictions even if
provided with a winner-takes-all operation. a, Outcomes of PredPCA offer long-term prediction via greedy prediction based on iterative winner-takes-all
operations, regardless of training dataset. Each row indicates a prediction based on a different realization of training sequence. A transition mappingfrom xt|t-1 to xt+1|t is assumed. b, The long-term prediction is successful even if a transition mapping from xt|t-1 ⊗ xt-1|t-2 to xt+1|t is assumed, indicating theminimal influence of the assumed model structure (that is, prior knowledge). c, PredPCA can also predict Fibonacci sequences in the long term, regardlessof the training dataset. d, Model selection to determine the optimal number of step backs. Here, the standard AIC was used for model selection. We
considered the following four models based on four types of polynomial basis functions, xt|t-1, xt|t-1 ⊗ xt-1|t-2, xt|t-1 ⊗ xt-1|t-2 ⊗ xt-2|t-3, and xt|t-1 ⊗ xt-1|t-2 ⊗ xt-2|t-3⊗ xt-3|t-4. The state in the next time period xt+1|t was predicted based on these four types of bases, followed by a winner-takes-all operation to conduct thegreedy prediction, and their AICs were compared. Left panel: To explain the ascending order sequences, a mapping from xt|t-1 to xt+1|t was the best among
these four models. Right panel: To explain the Fibonacci sequences, a mapping from xt|t-1 ⊗ xt-1|t-2 to xt+1|t was significantly better than other three models.
Here, the pairwise t test was applied based on 10 different realizations. Error bars indicate the standard deviation. e, In contrast, SSM based on Kalmanfilter tends to fail iterative prediction depending on the initial conditions of state and parameter values, and training history—even though it uses thewinner-takes-all operation—owing to its relatively large state and parameter estimation errors. System identification using SSM is severely harmed by
nonlinear interaction between state and parameter estimations, which yield local minima or spurious solutions (Extended Data Fig. 2b); consequently,
SSM exhibits an approximately 6% categorization error (Fig. 2b). These inaccuracies undermine iterative predictions using SSM even when states are
de-noised in each step using a winner-takes-all operation.

Nature Machine Intelligence | www.nature.com/natmachintell

Articles

NATurE MACHInE InTEllIGEnCE

Extended Data Fig. 4 | Instability of features extracted by TAE and SSM. This figure is a counterpart of Fig. 3b. TAE and SSM do not guarantee the global
convergence of their outcomes, and as a result their extracted features are sensitive to the initial conditions, order of supplying mini batches, and levelof observation noise. The extracted features in six trials are shown; the last three are outcomes trained with a large noise. The same training dataset wasused for all trials. However, as initial parameter values for TAE and SSM and order of supplying mini batches were varied, different features were extracted.
The difference in the observation noise level also altered their outcomes. These results imply the unreliability of features extracted by TAE and SSM, andfurther highlight the benefit of the global convergence guarantee of PredPCA.

Nature Machine Intelligence | www.nature.com/natmachintell

NATurE MACHInE InTEllIGEnCE

Extended Data Fig. 5 | See next page for caption.

Nature Machine Intelligence | www.nature.com/natmachintell

ArticlesArticles

NATurE MACHInE InTEllIGEnCE

Extended Data Fig. 5 | Feature extraction of diving car movies. a, PC1–PC3 of the categorical features (that is, x̄t) representing the brightness and verticaland lateral symmetries of scenes. b, PC1 of the dynamical features (that is, Δxt+3|t) representing the lateral motion. Although (a)(b) were obtained using
PredPCA with grouping of the data, these extracted features accurately matched those obtained using PredPCA without the six sub-groups (Fig. 4b,c).
This implies that PredPCA offers reliable identification of relevant features, even when using the data grouping. c, 100 major categorical features ( x̄t)
representing different categories of scenes. d, 100 major dynamical features (Δxt+3|t) responding to motions at different positions of the screen. The whiteareas indicate the receptive field of each encoder. c and d were obtained using PredPCA and ICA without the six sub-groups. Similar to Fig. 3b, theseimages visualize linear mappings from each independent component to the observation.

Nature Machine Intelligence | www.nature.com/natmachintell