Neural Networks 25 (2012) 191–199

Contents lists available at SciVerse ScienceDirect

Neural Networks
journal homepage: www.elsevier.com/locate/neunet

Real-time stylistic prediction for whole-body human motions
Takamitsu Matsubara a,b,∗ , Sang-Ho Hyon b,c , Jun Morimoto baGraduate School of Information Science, NAIST, 8916-5, Takayama-cho, Ikoma, Nara, 630-0101, Japanb

Department of Brain Robot Interface, ATR-CNS, 2-2-2, Hikaridai, Seika-cho, Soraku-gun, Kyoto, 619-0288, Japanc

Department of Robotics, Ritsumeikan University, 1-1-1, Nojihigashi, Kusatsu, Shiga, 525-8577, Japan

articleinfo

Article history:
Received 13 April 2011
Received in revised form 7 July 2011
Accepted 29 August 2011
Keywords:
Human motion prediction
Multifactor state-space models
Style-content separation
Stylistic prediction

abstract
The ability to predict human motion is crucial in several contexts such as human tracking by computer
vision and the synthesis of human-like computer graphics. Previous work has focused on off-line
processes with well-segmented data; however, many applications such as robotics require real-time
control with efficient computation. In this paper, we propose a novel approach called real-time stylistic
prediction for whole-body human motions to satisfy these requirements. This approach uses a novel
generative model to represent a whole-body human motion including rhythmic motion (e.g., walking)and discrete motion (e.g., jumping). The generative model is composed of a low-dimensional state (phase)
dynamics and a two-factor observation model, allowing it to capture the diversity of motion styles
in humans. A real-time adaptation algorithm was derived to estimate both state variables and styleparameter of the model from non-stationary unlabeled sequential observations. Moreover, with a simple
modification, the algorithm allows real-time adaptation even from incomplete (partial) observations.
Based on the estimated state and style, a future motion sequence can be accurately predicted. In our
implementation, it takes less than 15ms for both adaptation and prediction at each observation. Our realtime stylistic prediction was evaluated for human walking, running, and jumping behaviors.
© 2011 Elsevier Ltd. All rights reserved.

1. Introduction
Over the last decade, a considerable number of studies havebeen conducted on learning the generative models of humanmotion for modeling, prediction, and recognition (Howe, Leventon,& Freeman, 2000; Li, Wang, & Shum, 2002; Ormoneit, Sidenbladh,
Blank & Hastie, 2001; Pavlovic, Rehg & MacCormick, 2000;
Sidenbladh, Black & Fleet, 2000; Urtasun, Fleet, & Fua, 2006;
Urtasun, Fleet, Hertzmann & Fua, 2005; Wang, Fleet, & Hertzmann,2008). A significant limitation of these methodologies is that theycannot explicitly consider the natural variations of human motionsin the generative model, widely referred to as style (Brand &Hertzmann, 2000; Grochow, Martin, Hertzmann & Popovic, 2004;
Hsu, Pulli, & Popovic, 2005; Shapiro, Cao, & Faloutsos, 2006; Taylor& Hinton, 2009; Torresani, Hackney & Bregler, 2006; Wang, Fleet, &Hertzmann, 2007). For example, as illustrated in Fig. 1, even for anindividual, each walking motion sequence has a distinct walking
style. These differences can be much larger between different
individuals. Therefore, to achieve highly accurate prediction for
a newly observed motion sequence, adaptation of the generative∗ Corresponding author at: Graduate School of Information Science, NAIST, 89165, Takayama-cho, Ikoma, Nara, 630-0101, Japan.
E-mail address: takam-m@is.naist.jp (T. Matsubara).
0893-6080/$ – see front matter © 2011 Elsevier Ltd. All rights reserved.
doi:10.1016/j.neunet.2011.08.008

model to the motion sequence by capturing the style of thesequence is necessary.
While most previous studies have focused on off-line processes
with well-segmented data, many robotics applications (e.g.,human–robot interaction (Onishi, Luo, Odashima, Hirano, Tahara& Mukai, 2007), imitation learning by humanoids (Ijspeert,
Nakanishi, & Schaal, 2002; Inamura, Toshima, & Nakamura, 2002;
Riley, Ude, Wada, & Atkeson, 2003) and powered suits (Fukuda,
Tsuji, Kaneko, & Otsuka, 2003; Kawamoto, Kanbe, & Sankai,
2003)) require real-time control with high accuracy and efficientcomputation in the prediction procedure.
In this paper, we propose a novel approach called real-time
stylistic prediction for whole-body human motions. Unlike previous
studies (Brand & Hertzmann, 2000; Taylor & Hinton, 2009; Wanget al., 2007), as illustrated in Fig. 2, in our approach the generativemodel adapts to a newly observed motion sequence by estimating
its style by a real-time process. Being able to perform this process
in real-time is based on (1) the simple structure of the generativemodel and (2) the adaptation algorithm which requires small
computational effort. We propose a generative model for wholebody human motion that is composed of a low-dimensional state(phase) dynamics and a two-factor (phase dependent observation
bases and style parameter) observation model to capture thediversity of motion styles in humans. We also present a learning
procedure to acquire the model from a variety of motion sequences

192

T. Matsubara et al. / Neural Networks 25 (2012) 191–199

Fig. 1. Illustration of style in human motion sequences. Ten walking phase-aligned
sequences by two individuals are overwritten in order of phase. The style in walking
behavior is considered as a control variable for the spatial variations.

including a diversity of motion styles. A real-time adaptation
algorithm was derived using an on-line Expectation-Maximization(EM) algorithm for computationally efficient inference of both thecorresponding state variables and the style parameter from nonstationary unlabeled sequential observations. Moreover, with asimple modification, the algorithm allows real-time adaptationeven from incomplete (partial) observations. Such applicability of
the adaptation algorithm for partial observations is very importantin a practical sense because we often meet situations wheresome elements of the observations are missing due to the limited
number of sensors available or occlusions (Chai & Hodgins, 2005).
Based on the estimated state and style, the generative model canaccurately predict a future motion sequence.
On the other hand, most of the existing models that explicitly
estimate the style of motion can achieve neither real-time adaptation nor non-stationary motion estimation since the inference algorithm requires large computational effort (Brand & Hertzmann,
2000; Ormoneit et al., 2001; Sidenbladh et al., 2000; Taylor & Hinton, 2009; Urtasun & Fua, 2004; Wang et al., 2007).
The organization of this paper is as follows. In Section 2,
we present a novel generative model to represent a wholebody human motion including rhythmic motion (e.g., walking)and discrete motion (e.g., jumping). We also present a learning
procedure to acquire the model from a variety of motion sequences
including a diversity of motion styles. In Section 3, a real-time
adaptation algorithm is derived by applying an approximatedEM algorithm to the generative model. Moreover, we present a
simple modification in the adaptation algorithm, which allows
real-time adaptation even from incomplete (partial) observations.
A real-time prediction method of future motion sequence basedon the estimated state and style is also presented. In Section 4,
the effectiveness of our real-time stylistic prediction is validated
for human walking, running, and jumping behaviors with motion
capture data. Section 5 concludes this paper.
2. Learning generative modelThis section describes the proposed generative model and alearning procedure for the model with a stylistic data set.
2.1. Generative model for whole-body human motions
We first define the notation of the proposed generative model.
x ∈ Rd is the state variable, y ∈ RD is the observation and the probability distribution p(yt |xt ; w) is the observation model. p(xt +1 |xt )is the state-transition probability distribution. The parameter vector w ∈ RJ is an additional latent variable that controls the spatial variation of observations. We call this the style parameter. Itsgraphical model is depicted in Fig. 3. For periodic and discrete motions, we explicitly define the state variable x as
xt = [φt ωt ] =T[ψt ψ̇t ]T[pt ṗt ]T(Rhythmic)(Discrete).

observation
Estimation

Walking

Running

Striding

Real-time Adaptation for Style non-Stationary Motion SequenceFig. 2. Illustration of the real-time adaptation and prediction of the generativemodel for a non-stationary motion sequence with styles (walking behaviors). The
test sequence consists of three motions, walking, running and striding generated
by different individuals. The solid human figure is as observed and the dashed oneis the predicted motion as a result of adaptation of the generative model to theobservation sequence. The adaptation is achieved by on-line EM incrementally withlittle computation at each observation. For all motions, the model is rapidly adapted
to the style of the recent test sequence since the time-forgetting factor effectively
forgets past observations.

S ⊂ R2 and its velocity ω ≡ ψ̇ to represent its periodicity of rhythmic motions, similar to Ormoneit et al. (2001) and Urtasun and Fua(2004). We also define the phase φ as a point on a one dimensional
closed line segment φ ≡ p ∈ L for discrete motions to represent its
non-periodicity (discreteness). The explicit use of these assumptions in the generative model yields the low-dimensional statevariable x. Moreover, as presented in the next section, it allows a
simple learning algorithm for the generative model from data.
Based on the above assumptions, we conclude that the statetransition model and the observation model are modeled by
Gaussian distributions as:p(xt +1 |xt ) = N (µx (xt ), Σx (xt )),(2)p(yt |zt ; w) = N (µy (zt ; w), Σy (zt ; w)),

(3)

where
zt = g (xt ) =[cos(φt ) sin(φt )]Tφt(Rhythmic)(Discrete)(4)

and the observation model is defined as a probabilistic mappingfrom a phase φt (as zt ) to an observation yt . The velocity of phaseωt governs the temporal variation of the time-series, that is, it
controls the velocity of human motions generated by the model.
For rhythmic motions, zt = g (φt ) ∈ R2 represents a point on amanifold S in R2 where the radius is r = 1 and the angle is φt .
This state representation allows us to approximately measure thegeodesic distance between points on S as the Euclidean distance in
R2 . For discrete motions, zt is the equivalent of φt .1
2.2. Learning procedure with a stylistic data set
The learning procedure assumes we have multiple humanmotion sequences including a diversity of motion styles. Let Ys =[ys1 · · · ysC (s) ]T ∈ RC (s)×D denote a time-invariant motion sequencewith a distinct style, where s ∈ {1, 2, . . . , S } is the style index
in which each value indicates a corresponding distinct style, c ∈{1, 2, . . . , C (s)} is the content index that corresponds to the phase(1)

That is, we define the state variable x by phase φ as a point on a one
dimensional sphere in two dimensional Euclidean space φ ≡ ψ ∈1 For the case of a discrete motion, z is a scalar; however, it is kept as a vectortnotation for simplicity of the overall description.

T. Matsubara et al. / Neural Networks 25 (2012) 191–199

193

VT
Let Yallbe an Sa× DC matrix stacked from DS × C Matrix
Yall
a . Then, SVD for this matrix leads to the following factorial
representation as
Yall
a

VT

= USVT ≈ WỸ.

(6)
1

Fig. 3. A generative model representing time series data for a human motion. The
state variable is defined as xt = [φt ωt ] where the phase φ of the state variable
is defined as a point on one dimensional sphere in two dimensional Euclidean
space φ ≡ ψ ∈ S ⊂ R2 and as a point on one dimensional closed line segment
φ ≡ p ∈ L for rhythmic and discrete motions, respectively. The observation
yt is conditionally independent from all other variables given the state xt . xt +1
is conditionally independent given the state xt . w is the style parameter vector
invariant for the time.

and ysc

D∈ R is an observation with the style indexed by s andcontent c. In the following, let us assume that a set of training
sequences D = {Y1 , . . . , YS } is given for learning a generativemodel as a stylistic data set. The learning procedure is composed of
the following three steps: (1) Data alignment by phase information,(2) Extraction of observation bases, (3) Learning generative modelfrom the bases. These three steps achieve learning of a compact
generative model of a human motion that can represent large
stylistic variations, where both the state variables and the styleparameters of the generative model can be low-dimensional.
2.2.1. Data alignment by phase information
For the subsequent steps, we align all motion sequences in astylistic data set by phase. We introduce data alignment strategies
for rhythmic and discrete motions.
Rhythmic motions: We utilize auto-correlative and cross-correlative
coefficients for the alignment process. First, we maximize the autocorrelative coefficient for identifying the period T of each sequence∑N s T s
as: T s ← arg maxj As (j), where As (j) =n yn yn+j is the autocorrelative coefficient with the self-index shift j in phase with astyle indexed by s. Next, we maximize the cross-correlative coefficient to find the optimal cross-index shift h in phase as: hs ←arg maxj C s (j), where C s (j) =bT sn yn yj+rd( nT s ) , and index b is the∑N

Tbstyle index corresponding to the sequence that has the shortest period. T b is the period of the shortest period indexed by b. The function rd(·) is a round-off function. The above procedures yield analigned data matrix:

 1
y

1

rd(h1 + T b )



 S

Yall
a = 

y

..
.

T

S
rd(hS + T b )T···y1rd(h1 +T 1 )



..

..
.






.
···

ySrd(hS +T S )

(5)

S T
DS ×C
1
and Ysa ∈ RC ×D . Note that eachwhere Yalla = [Ya · · · Ya ] ∈ R
allrow of the matrix Ya are the observations corresponding to thesame value of the phase φ . Each column indicates a corresponding
motion sequence indexed by s.

Discrete motions: We align each motion sequence Ys from the samestarting point φ = 0 to the goal point φ = 1. The aligned datasmatrix Yalla can be obtained by simply setting T as the duration of
motion and hs = 0 in Eq. (5) for all s.
2.2.2. Extraction of observation bases
Since the aligned data matrix Yalla is a rectangular matrix, Singular Value Decomposition (SVD) based matrix factorization canbe applied to extract observation bases. By applying the factorization, we can form a style-content factorial model (referred to as theasymmetric bilinear model in Tenenbaum and Freeman (2000)).

S T

S ×J

We define the style parameter matrix W = [w · · · w ] ∈ R
to be the first J (≤S ) rows of U, and the content parameter matrix
Ỹ = ([Ȳ1 · · · ȲJ ]T )VT ∈ RJ ×DC to be the first J columns of SVT . As∑Js ja result, we can obtain an approximated form as Ysa ≈j=1 wj Ȳ .
The obtained matrix Ȳj ∈ RC ×D is named the j-th observation basis,and the vector ws ∈ RJ is the s-th style parameter vector.
2.2.3. Learning generative modelWith the extracted observation bases Ȳj for all j, we can learn
the generative model. Since each point of an observation basis
corresponds to a value of phase on S ⊂ R2 (or L), we learn amapping between ȳj and z using all data of each basis. Here weutilize Gaussian process regression (Rasmussen & Williams, 2006)
because it allows us to derive an analytically tractable predictive
distribution and to learn hyperparameters by maximization of themarginalized likelihood.
Each basis Ȳj is independently modeled as a Gaussian process
as





1
p(Ȳj |Z, βj ) ∝ exp − Tr((Kjy )−1 Ȳj (Ȳj )T )
2

(7)
j

where Z is the phase-aligned matrix corresponding to Ȳj and Ky ∈

R

C ×C

is the gram matrix in which (p, q)

j
1 exp

β

j

β2

j
entry is ky(zp , zq ) =j
3 , and the hyperparameter is(− 2 ‖zp − zq ‖ ) + δzp ,zq /βjjj
represented as βj = {β1 , β2 , β3 }. With the above settings, the2predictive distribution of the j-th observation basis ȳj∗ given a
novel input z∗ can be easily derived asp(ȳj∗ |z∗ , Ȳj , Z) = N (µ̄j (z∗ ), Σ̄ j (z∗ ))

(8)

where

µ̄j (z∗ ) = (Ȳj )T (Kjy )−1 kjy (z∗ ),(9)Σ̄ j (z∗ ) = (kjy (z∗ , z∗ ) − kjy (z∗ )T (Kjy )−1 kjy (z∗ ))I,j
and ky

j
ky

j
ky

(10)

(z ) = [ (z , z1 ) · · · (z , zC )] (Rasmussen & Williams,
2006). The predictive distribution for observation y∗ given a novel
input z∗ with a style parameter vector ws can be written withthe standard result of a Gaussian distribution and invariant linear
transformation as
∗

∗

∗

T

p(y∗ |z∗ , Ȳ1:J , Z; ws ) = N (µy (z∗ ; ws ), Σy (z∗ ; ws ))

(11)

where

µy (z∗ ; ws ) =J−

wjs µ̄j (z∗ ),

(12)

J
−
(wjs )2 Σ̄ j (z∗ ).

(13)

j =1

Σy (z∗ ; ws ) =j =1

This predictive distribution is set as the observation modelp(yt |xt ; w) of the generative model in Eq. (3).
The state transition model p(xt +1 |xt ) of the generative modelin Eq. (2) is simply modeled under the assumption that all data areobserved with a fixed sampling rate as p(xt +1 |xt ) = N (Axt , Q)
where the state transition matrix is A =

1
0

1
1

and the covariancematrix of process noise Q should be set appropriately for the databy hand. The mean and variance of the state transition model are
written as µx (xt ) = Axt , Σx (xt ) = Q, respectively.

194

T. Matsubara et al. / Neural Networks 25 (2012) 191–199

3. Real-time stylistic prediction

3.2. On-line implementation of the proposed estimation method

In this section, we derive a real-time adaptation algorithm forthe generative model to a newly observed test sequence by using
an approximated EM algorithm, i.e., the algorithm that estimatesa time-series of state variables X̂ = [x̂1 · · · x̂T ]T and the styleparameter vector ŵ from a portion of the test sequence Ŷ =[ŷ1 · · · ŷT ]T in an on-line manner. In Section 3.1, we start from thederivation of off-line (batch) computations of E-step and M-stepby introducing simple approximations. This gives us an efficientcomputation for the adaptation of the generative model for atest sequence. In Section 3.2, the derived algorithm is further
extended to a fully recursive, time-forgetting and computationallyefficient on-line algorithm by a modification of the likelihood assuggested by Sato and Ishii (2000). Section 3.3 describes how we
can achieve real-time stylistic prediction of a test sequence withthe generative model and the adaptation algorithm. Section 3.4
presents a modification of the adaptation algorithm to allow realtime adaptation even from incomplete (partial) observations.

Here we present an on-line implementation of the proposed
estimation method. As suggested by Sato and Ishii (2000), we
introduce a time-forgetting factor in our likelihood in Eq. (15).
This yields a recursive, time-forgetting on-line EM algorithm, inwhich both the E-step and M-step can be calculated at each newobservation with little computation.
In the on-line EM algorithm, the weighted mean is defined as

3.1. Estimation of state variables and style parameters by using anapproximated EM algorithm

where ηT = {1 + λT /ηT −1 } . Introducing the above weighted
mean in the computation of the EM algorithm yields an on-line
EM algorithm. i.e., a recursive, time-forgetting on-line adaptation
algorithm. If we set λs = λ < 1.0, this can be applied for nonstationary motion sequences with styles. The details of the derived
on-line adaptation algorithm are presented in Appendix.

To derive an efficient adaptation algorithm for the generativemodel to a test sequence, we fit the EM algorithm to the
problem with simple approximations. The EM alternates between
optimizing a distribution over the state variables (the E-step) and
optimizing the style parameters given the posterior distribution
obtained by the E-step (the M-step) (Dempster, Laird, & Rubin,
1977; Ghahramani & Hinton, 1996; Shumway & Stoffer, 1982).
Introducing a distribution over state variables as Q(X̂), a lower
bound of the log-likelihood F (Q, ŵ) given observations Ŷ can bederived from the likelihood L(Ŷ|ŵ). The E-step holds the styleparameters fixed and sets Q to be the posterior distribution overthe state variables:

Qk+1 (X̂) ← arg max F (Q, ŵk ) = p(X̂|Ŷ, ŵk ).

(14)

Q

This maximizes F w.r.t. Q turning the lower bound of the
log-likelihood into the likelihood as F (Q, ŵ) = L(Ŷ|ŵ). The Mstep holds the distribution fixed and computes the style parameters that maximizes F as ŵk+1 ← arg maxŵ F (Qk+1 , ŵ) =
arg maxŵ X̂ Qk+1 (X̂) log p(X̂, Ŷ|ŵ)dX̂. Therefore, the E-step ink-th iteration calculates the posterior of the state variable X asp(X̂|Ŷ, ŵk ), and the M-step updates the new style parameter ŵk+1
with the posterior. This alternative procedure typically converges
to a locally optimal point through a small number of iterations.
Analytical computations of the E-step and M-step for our generative model, however, are intractable due to the nonlinearity inobservation bases against phase. Here, we introduce simple approximations to create analytical solutions for both steps. Under
the assumption that all observation bases are smooth functions
over the phase variable φ , the E-step can be approximately calculated in the same way as the Extended Kalman Filter (EKF) (Ko &Fox, 2008). Thus, we approximate Q ≈ N (X̂) to obtain an analytical computation of the E-step. Subsequently, the analytical computation of the M-step is derived with the following approximation as∫p(X̂|Ŷ; ŵk ) log p(X̂, Ŷ; ŵ)dX̂

ŵk+1 ← arg max
ŵ≈ arg max
ŵ

X̂
T
−

log p(ŷt |x̂pm,t ; ŵ)

(15)

t =2

where X̂pm is the mean of the posterior distribution as X̂pm =X̂p(X̂|Ŷ; ŵk )dX̂ which is obtained by the E-step. The details of
the derived adaptation algorithm are presented in Appendix.

X̂

⟨⟨f (x)⟩⟩T ≡ ηT

T
−
t =1



T
∏


λs f (xt )(16)

s=t +1

where ηT ≡ { t =1 { s=t +1 λs }}−1 and the parameter λs (0 ≤ λs≤ 1) is a time-dependent forgetting factor which is introduced forforgetting the effect of the past observations. ηT is a normalized
coefficient and plays a role similar to a learning rate. This weighted
mean ⟨⟨·⟩⟩i has a step-wise equation∑T∏T⟨⟨f (x)⟩⟩T = (1 − ηT )⟨⟨f (x)⟩⟩T −1 + ηT f (xT )

(17)

−1

3.3. Stylistic prediction of future observations
With the generative model in which both the state variables andthe style parameter are estimated from a given test sequence, we
could predict its future states and observations by the following
simple algorithm. Given a portion of a test sequence until time t, we
can estimate the expectation of the corresponding state variablex̂t and the estimated value of the style parameter ŵt . One-stepahead prediction of the expectation of state variable x̂t +1|t can beobtained by x̂t +1|t = p(xt +1 |x̂t )xt +1 dxt +1 = Ax̂t . Corresponding
estimation of observation ŷt +1|t is then obtained from ẑt +1|t =g (φ̂t +1|t ) and Eq. (3) as ŷt +1|t = p(yt +1 |ẑt +1|t ; ŵt )yt +1 dyt +1 =µ(ẑt +1|t ; ŵt ). With these calculations, stylistic prediction byidentifying the style of the test sequence can be achieved.
Executing the prediction recursively, it is possible to achieve
multiple-step ahead predictions of future states and observations.



3.4. Adaptation from incomplete (partial) observations
In this section, we present a modification in the adaptation
algorithm so that it allows real-time adaptation even fromincomplete (partial) observations. In a practical sense, suchapplicability of the adaptation algorithm to partial observations
is very important because we often meet situations where someelements of the observations are missing due to the limited
number of sensors available or occlusions. Our approach is inspiredby a model for dealing with packet losses in a wireless networkwith the Kalman filter (Liu & Goldsmith, 2004).
For simplicity (but without loss of generality), we divide anobservation into two parts as y = [y̌T , ỳT ]T where y̌ ∈ RD−Dm areobservable elements and ỳ ∈ RDm are missing observations, that is,
we assume that the index of missing elements in the observation
is known. To manage such missing observations for the adaptation
algorithm, we modify the observation model in Eq. (3) asp

[ ] 

y̌ 
z
;
w
ỳ  tt[] [µ̌y (zt ; w)
R
=N
, 11
µ̀y (zt ; w)
R21

R12

R22 + σp Im

]
(18)

T. Matsubara et al. / Neural Networks 25 (2012) 191–199

where µy = [µ̌y , µ̀y ]T is the mean, R =T

T(D−Dm )×(D−Dm )



R11
R21

R12
R22



is the co-

RT21 is R(D−Dm )×Dm

variance matrix, R11 is R
, R12 =
and R22 is RDm ×Dm matrices. As indicated in Liu and Goldsmith(2004), in the above model, the missing observations ỳ can be represented by setting σp → ∞ in Eq. (18). With this setting, an
adaptation algorithm for the partial observation can be derived byintroducing the following modifications in the fully observable
case (see Appendix) as
y → y̌,µy → µ̌y ,(19)

H → Ȟ,R → R11 (20)

where H =ȞH̀and Ȟ is R(D−Dm )×Dm , H̀ is RDm ×Dm matrices. The

derivation is based on a formula of the inverse of a partitioned matrix (Liu & Goldsmith, 2004).
The derived adaptation algorithm allows us to adapt thegenerative model of whole-body human motions to a test sequence
from partially observed information.
4. Experimental results
In this section, we validate the effectiveness of our realtime stylistic prediction method for human motions using data
captured by a motion capture system. Section 4.1 first presents theeffectiveness for rhythmic motions such as walking and runningbehaviors in comparison with several on-line prediction methods
with respect to computational costs and prediction accuracy. Its
real-time calculability and applicability for partial observations arealso validated through experiments. In Section 4.2, the proposed
method is applied to discrete motions such as forward andvertical jumping behaviors to demonstrate the effectiveness of the
proposed method for several kinds of behaviors.
4.1. Stylistic prediction for walking motions
We first learned a generative model from a collection of motion sequences. The data were taken from the CMU Graphics Lab’s
motion capture database (http://mocap.cs.cmu.edu). Each observation yt = [qTt rTt vTt ]T ∈ R62 was down-sampled to 60 Hz, andfocused on full-body joint Euler angles qt ∈ R56 . To obtain a diversity of walking styles for this experiment, we selected several
motion sequences of walking, running and striding motions captured from four subjects (subject-A, B, C, D).2 More specifically,
we selected eight walking sequences and one striding sequence of
subject-A, four running sequences of subject-B, two running sequences and one walking sequence of subject-C, and one runningsequence of subject-D.3 The training sequences for model learning
were chosen as six walking sequences and one striding sequenceof subject-A, three running sequences of subject-B and two running sequences of subject-C. The dimension of the style parameter
vector was set as J = 4 with a guide of the spectrum of singularvalues obtained. The estimated observation bases with the above
settings are plotted in Fig. 4(e).
4.1.1. Verification of the learned generative modelFirst, we verified the learned generative model for stylistic
prediction using the batch EM algorithm. For several segmented

2 The correspondence of the subjects to the labels in the CMU motion capture
database is subject-A:08 subject-B:35, subject-C:02 and subject-D:16.
3 {08_01–08_09, 35_20–35_23, 02_01–02_03, 16_35}.amc. Training sequences
contain {08_02–08_08, 35_20–35_22, 02_02–02_02 }.amc.

195

test sequences, the prediction accuracy of the adapted generativemodel was evaluated. The adaptation was executed with the EM
algorithm using a portion of the test sequence, which was set as
0.0–2.0 s (120 frames) for walking and 0.0–1.0 s (60 frames) for
running, based on their motion frequencies. Another portion for
0.5 s (30 frames) was then used as the ground truth for evaluationof the prediction accuracy. Test sequences were set as a striding
sequence of subject-A, a running sequence of subject-B, a walking
sequence of subject-C, and a running sequence of subject-D. Notethat none of the test sequences were included in the training
sequences, and no motion sequence of subject-D was included in
the training sequences, i.e., subject-D was a completely unknown
subject for the generative model. As a criterion for evaluation of
the predictive accuracy, we used the average prediction error:
∑G+i
1
off
2
Erms
= GD
t =i+1 ‖yt − ŷt |i ‖ , where G is the number of predicted
frames, D is the dimension of observation and i is the index
of the terminal point of the test sequence used for adaptation.
In all cases, the EM procedure converged within 10 iterations.
Fig. 4(a)–(d) illustrates the result of the stylistic prediction for thetest sequences compared with the ground truth. Each predicted
pose is very similar to the ground truth for all cases including
the unknown subject. The estimated style parameters are plottedin Fig. 4(f). As shown in (f), the estimated style parameters werelargely different due to the distinct style of each test sequence,
producing high accuracy predictions. The average prediction error
off
for all test sequences Erms
= 1.68 with J = 1 is 58% reduced byoffusing J = 4 as Erms = 0.71. These results verified the learned
generative model for stylistic prediction.
4.1.2. Real-time stylistic prediction
Next, we evaluated our real-time stylistic prediction methodwith the derived adaptation algorithm and the learned generativemodel. As a non-stationary motion sequence with styles, a test
sequence was prepared successively uniting three different test
sequences. In this experiment, the adaptation procedure andprediction at each observation approximately required less than
15 ms in our Matlab implementation on a PC (Intel(R) Core-i7 CPU
3.33 GHz), i.e., about 13 ms for the E-step, 0.2 ms for the M-stepand 0.5 ms for one-step ahead prediction of observation. Comparedwith the average walking period T = 1.2 s this processing
time would be relatively short and effective enough for real-time
prediction.
We compared the prediction performance of our proposed
method with standard on-line prediction methods. Specifically, we
compared it with recursive least square linear regression (RLS)
with a time-forgetting factor λ (Haykin, 2002) and with Radial
Basis Function Networks (Bishop, 1995) with a time-forgetting
factor as standard approaches for on-line function approximation.
For the RBFs, principal component analysis (PCA) was applied
separately for each test sequence and obtained a low-dimensional
subspace (RBFs + PCA). The basis was located on a grid with aneven interval in each dimension of 3-dimensional feature spaceobtained by PCA. We used 1000 (=10 × 10 × 10) basis functions inthis experiment. As a nonparametric approach, we implemented aparticle filter as used in Sidenbladh et al. (2000) and Ormoneit et al.
(2001). In this approach, the required posterior density function
is represented by a set of random samples {[xk wk ]T }Kk=1 withkTassociated weights mk as p(xt , wt |y1:T ) ≈k=1 m δ([x w] −kk T
[xt wt ] ) where K is the number of particles. As the numberof samples becomes very large, this approximation becomes anequivalent representation of the posterior in the sense of an
optimal Bayesian estimate. The proposal distribution for random
sampling of particles and the likelihood model were set by thestate-transition model and the observation model of our proposedmodel as in Eq. (3). The transition model of the style parameter was
newly set as wt +1 ∼ N (wt , σw I).

∑K

196

T. Matsubara et al. / Neural Networks 25 (2012) 191–199

e

1

Y

a

Phase Dependent Basis

b

2

Y

3

Y

4

Y

0

2π

φ

f

d

Style Parameters

c

π
Phase

w1
w2
w3
w4

0.75
0.5
0.25
0.0
-0.25

Running
Striding
Walking Running(Subject-A) (Subject-B) (Subject-C) (Subject-D)

Fig. 4. Results of off-line stylistic prediction for human walking and running sequences. (a)–(b) show the prediction results where solid-lines indicate the ground truth ineach frame and dash-lines show the predicted pose. The observation bases are shown in (e). The estimated style parameter w for test sequences is plotted in (f) where each
wj corresponds to Ȳj for all j. (a) Shows the prediction results for a striding-style walking sequence, (b) shows the results for a running sequence and (c) shows a normal
walking sequence. These sequences seem to have different motion styles and such differences in motion style are captured by the estimated style parameters as shown in(f) to achieve the predictions with high accuracy. (d) Also shows the results for a running sequence which seems to be a slightly different motion style from (b). Such a small
difference in motion style is also captured by the estimated style parameters as in shown (f).
Table 1
Prediction results for non-stationary human walking data in our Matlab implementation on a PC (Intel(R) Core-i7 CPU 3.33 GHz).
MethodTime (s)on
Erms

Proposed
RLS
RBFs + PCA
PF (K = 500)PF (K = 10,000)

0.01
2.00 × 10−4
0.09
0.13
105

0.85
1.76
1.40
1.73 (0.05)
1.23 (0.05)

The performance of the stylistic prediction was evaluated by∑∑ G +i1onthe on-line average prediction error Erms
= T1 Ti=1 { GDt =i+1 ‖yt − ŷt |i ‖2 }, where G is the number of predicted frames, D isthe dimension of observation, T is the number of total frames of
the test sequence and i is the index of the observed frame. The
predicted pose ŷi|i is ŷi|i = µ(x̂i ; ŵi ) where, x̂i and ŵi are estimatedby the adaptation algorithm using sequential observation until yi .
Table 1 shows the average prediction error in predictions (G =25 frames). At each observation for a test sequence, adaptationand prediction were executed. Then, the prediction performanceon
was evaluated by the average prediction error Erms
for all cases.
For the proposed method, the time-forgetting factor was set as
λ = 0.8 so that the learning algorithm can quickly adapt to newly
acquired data while old data is properly discarded. For the RLS
and RBFs + PCA, the forgetting factor was set as λ = 0.95. Since
the number of parameters required in these methods are muchlarger than the proposed method, the larger forgetting factor isnecessary to properly evaluate the likelihood for the parameters.
From Table 1, the proposed method showed the best performance
among all methods in terms of both computational cost andprediction accuracy. The results of these comparisons suggestthat it is impossible to represent human motion by well-known
linear dynamics. Moreover, while the use of non-linear dynamics
represented by RBFs potentially captures the non-linearity; thehigh-dimensionality resulting from the use of non-linear function
approximators makes it difficult to achieve successful modelingwith a modest amount of data obtained in on-line. A particle
filtering approach can be competitive in terms of prediction

Table 2
Prediction from partial observations.
Missing portions

on
Erms

(1) One arm(2) One leg(3) One leg and one arm

0.87
0.91
0.91

performance with a huge number of particles; however, it requires
prohibitive computational cost for real-time processes.
4.1.3. Prediction from partial observation
We also evaluated our prediction method for partial observations. For this experiment, three conditions were considered(1) missing one arm, (2) missing one leg, and (3) missing both one
arm and one leg.4 For such scenarios, the proposed method was applied for predictions and the results evaluated by the on-line prediction error using the non-stationary motion sequence with styles
as a test sequence. Table 2 shows the experimental results. In allcases, the prediction performance was not much impaired by missing data. These results suggest the effectiveness of the proposed
method for partial observations. Fig. 5 depicts prediction resultsfor partial observations.
For further validation of the proposed method for partial
observations, we also applied the method to a different data
set from the CMU motion capture database, which has beenused in Lawrence (2007) and Taylor, Hinton, and Roweis (2006).
We prepared 31 motion sequences of one subject (composed ofwalking and running) as training data and also prepared one
walking sequence as test data.5 As in Lawrence (2007), the test
sequence was modified in the two ways: (1) missing right leg and(2) missing upper body. The prediction accuracy was evaluated by

4 One arm includes 7 DoFs: radius(1), wrist(2), hand(1), finger(1), thumb(2). One
leg includes 7 DoFs: femur(3), tibia(1), foot(2), toes(1).
5 The training data was composed of {35_1–17, 19–26 and 30–34}.amc, and the
motion sequence 35_29.amc was used as test data.

T. Matsubara et al. / Neural Networks 25 (2012) 191–199

197

Root Mean Square Error

14aRight Leg
Upper Body

12
10
8
6
4
0

b

0.1

0.2

0.3

0.4
0.5
Time [s]

0.6

0.7

0.8

Fig. 6. Results of the on-line predictions of missing body movements. The graphindicates the root mean square error of the missing joint angles averaged over 60
future frames at each observation. ‘‘right leg’’ and ‘‘upper body’’ indicate the thecases of (1) missing right leg and (2) missing upper body in the test sequence,
respectively. For the both cases, our method rapidly converged and resulted in goodprediction accuracy of missing body movements.

aFig. 5. Results of predictions for walking behaviors from partial observations (Leftleg information is missing). Solid-lines indicate the ground truth in each frameand dash-lines show the predicted pose. (a) Shows the result for a striding-style
walking sequence and (b) shows the result for a running sequence. For both cases,
even though the style of each sequence is largely different and left-leg information
is missing, the predicted poses are very similar to the ground truth.

root mean square error of the missing joint angles averaged over
60 frames.
Fig. 6 shows the evaluation results of the on-line predictions.
Each line was obtained by applying a zero-phase shift low-passfilter with 1.5 cut-off frequency to the result of the on-line
predictions. The dimension of the style parameter vector was set
as J = 2, and the time-forgetting factor was set as λ = 0.8. Fig. 6
indicates that, for both (1) and (2), our method rapidly convergedat around 0.4 s and resulted in good prediction accuracy of missing
body movements. As presented in Lawrence (2007) and Taylor
et al. (2006), for this data set, the prediction of the missing upper
body is easier than that of the missing lower limb. Note that only
our proposed method is able to achieve the adaptation and theprediction on-line in real time.
4.2. Stylistic prediction for jumping motions
Next, we evaluated the proposed method for jumping motions.
Again, the data were taken from the CMU database with thesame setting of walking motions described in the previous section.
To include a diversity of jumping motion styles (e.g., forward
or vertical) in this experiment, we selected several motion
sequences of jumping motions captured from two subjects(subject-E, F). More specifically, we selected seven jumping
sequences of subject-E and two jumping sequences of subject-F.6
The dimension of the style parameter vector was set manually as
J = 6 by considering the range of singular values.
Using the learned generative model, we applied on-line
prediction to two test sequences. For validation, the first 60 frames
were utilized for on-line adaptation of the style and state variables,

6 The correspondence of the subjects to the label in the CMU motion
capture database is subject-E:13 subject-F:16. Training sequences contain{13_11, 13_13, 13_19, 13_39–13_42, 16_01–16_02}.amc. Test sequences contain{13_32, 13_40}.amc.

bFig. 7. Results of predictions for jumping motions. Two test sequences were usedfor validation. Initial motion sequences before takeoff are used for adaptation, andthen future motion sequences after takeoff are predicted and compared with theground truth to calculate prediction errors. Solid-lines indicate the ground truthin each frame and dash-lines show the predicted pose. (a) Shows the result fora vertical jumping motion after a deep squat motion. (b) Shows the result for a
forward jumping motion with a large arm swing. The proposed method captured
these motion-sequence specific features well in the predicted poses.

and then future observations were predicted and validated based
offon the average prediction errors Erms
(G = 30 frames) for bothcases. As a result, the errors were 0.97 and 1.05, and these values
were significantly small compared with the standard deviationsof both test sequences at 6.09 and 3.50. Fig. 7 depicts snapshotsof both predictions and the ground truth. For both cases, eachpredicted pose is very similar to the ground truth. More concretely,
in Fig. 7, (a) shows the result for a vertical jumping motion where
the hips and knees are used for a vertical jump, while (b) shows theresult for a forward jumping motion where the arms swing morethan the hips and knees. The proposed method captured thesemotion-sequence specific features very well in the predicted poses.
These results verified the effectiveness of our prediction method
even for discrete motions, and suggests the applicability of thismethod for a wide range of human motions.

198

T. Matsubara et al. / Neural Networks 25 (2012) 191–199

5. Conclusions

M-stepIn this paper, we proposed a novel approach called real-time
stylistic prediction for modeling and predicting human motions.
Our algorithm can be run in real-time and manage non-stationary
motion sequences with styles by capturing the style on-line withrelatively little computation. Moreover, the algorithm allows realtime adaptation even from incomplete (partial) observations. Its
effectiveness was demonstrated on motion capture data.
In the proposed approach, human motions are modeled by agenerative model which is composed of a low-dimensional state(phase) dynamics and a two-factor (phase dependent observation
bases and style parameter) observation model designed for
capturing the diversity of motion styles in humans. To achieve
higher prediction accuracy, a more complex model as the state
dynamics (e.g., nonlinear state transitions or factorized by multiple
variables) would be necessary. The proposed approach can
incorporate such a complex state dynamics by extending both theadaptation and prediction algorithms. For example, a nonlinear
state transition model can be incorporated in the proposed
approach. In such a case, the adaptation algorithm is derived withthe linearization of the state transition about the current meanof the state variable, and the prediction algorithm is executed by
nonlinear regression at each time evolution. Thus, with complex
state dynamics, the adaptation and prediction algorithms require
additional computational effort which makes it difficult to be
implemented in real time. Therefore, there is the trade-off between
the prediction accuracy and the computational effort required in
adaptation and prediction algorithms.
The ability to predict human motion is crucial in several
contexts such as human tracking by computer vision and thesynthesis of human-like computer graphics. Previous work has
focused on off-line processes with well-segmented data; thus,
many applications such as robotics require real-time control withefficient computation. The proposed method could be applied
to these applications. Our future work will address real robotics
applications such as human–robot interaction, imitation learningby humanoids and powered suits.
Our algorithm has free parameters such as time-forgetting anddimension of style parameter. Our future work will also address
setting them by a relevant determination from data.

1 T
ŵt = ⟨vT v⟩−
t ⟨v y⟩t ,(A.7)

where,
vt = [µ̄1 (x̂t ) · · · µ̄J (x̂t )]

(A.8)

T

⟨·⟩T =

1−
T t =1(·t ).

(A.9)

On-line M-step−1ŵt = ⟨⟨µT µ⟩⟩t ⟨⟨µT y⟩⟩t

(A.10)T⟨⟨µ µ⟩⟩t = (1 − ηt )⟨⟨µ µ⟩⟩t −1 + ηt µ(x̂t ) µ(x̂t )(A.11)⟨⟨µ y⟩⟩t = (1 − ηt )⟨⟨µ y⟩⟩t −1 + ηt µ(x̂t ) yt


λt −1
ηt = 1 +η t −1(A.12)

T

T

T

T

T

(A.13)

whereµ(xt |t −1 ) = [µ1 (xt |t −1 ) · · · µJ (xt |t −1 )](A.14)µ(x̂t ) = [µ (x̂t ) · · · µ (x̂t )](A.15)

1

J

T⟨·⟩T =

1−(·t )T t =1⟨⟨·⟩⟩T = ηT

T
−
t =1


ηT =

T
−
t =1





(A.16)T∏λs (·t )(A.17)

s=t +1

T
∏

−1
λs

.

(A.18)

s=t +1

Appendix B. Supplementary data
Supplementary material related to this article can be foundonline at doi:10.1016/j.neunet.2011.08.008.
References

Acknowledgments
This research was supported by the Strategic Research Programfor Brain Sciences ‘‘Brain Machine Interface Development’’ by the
Ministry of Education, Culture, Sports, Science and Technology,
Japan, and also partially supported by the Grant-in-Aid for
Scientific Research from Japan Society for the Promotion of Science(WAKATE-B22700177).
Appendix A. Derived EM steps
E-stepxt |t −1 = Ax̂t −1(A.1)

Σx,t |t −1 = AΣ̂x,t −1 AT + Q∂µ(x; ŵt −1 ) 
Ht |t −1 =∂xx=xt |t −1(A.2)

Kt = Σx,t |t −1 HTt|t −1 · (Ht |t −1 Σx,t |t −1 HTt|t −1 + R)

(A.4)

x̂t = xt |t −1 + Kt (yt − µ(xt |t −1 ; ŵt −1 ))(A.5)Σ̂x,t = (I − Kt Ht |t −1 )Σx,t |t −1 .

(A.6)

(A.3)

Bishop, C. M. (1995). Neural networks for pattern recognition. Oxford University Press.
Brand, M., & Hertzmann, A. (2000). Style machines. In SIGGRAPH (pp. 183–192).
Chai, J., & Hodgins, J. K. (2005). Performance animation from low-dimensional
control signals. ACM Transactions on Graphics, 24, 686–696.
Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood fromincomplete data via the EM algorithm. Journal of the Royal Statistical Society,
Series B, 39, 1–38.
Fukuda, O., Tsuji, T., Kaneko, M., & Otsuka, A. (2003). A human-assisting manipulator
teleoperated by EMG signals and arm motions. IEEE Transaction on Robotics andAutomation, 19, 210–222.
Ghahramani, Z., & Hinton, G. E. (1996). Parameter estimation for linear dynamicalsystems. Technical report. Dept. Computer Science, Univ. of Toronto (pp. 1–6).
Grochow, K., Martin, S. L., Hertzmann, A., & Popovic, Z. (2004). Style-based inverse
kinematics. ACM Transactions on Graphics, 23(1), 522–531.
Haykin, S. (2002). Adaptive filter theory. Prentice Hall.
Howe, N. R., Leventon, M. E., & Freeman, W. T. (2000). Bayesian reconstruction of
3D human motion from single-camera video. In Advances in neural information
processing systems: Vol. 12 (pp. 820–826).
Hsu, E., Pulli, K., & Popovic, J. (2005). Style translation for human motion. ACM
Transactions on Graphics, 24(3), 1082–1089.
Ijspeert, A. J., Nakanishi, J., & Schaal, S. (2002). Learning attractor landscapes forlearning motor primitives. In Advances in neural information processing systems:
Vol. 15 (pp. 1523–1530).
Inamura, T., Toshima, I., & Nakamura, Y. (2002). Acquisition and embodiment of
motion elements in closed mimesis loop. In IEEE international conference on
robotics and automation (pp. 1539–1544).
Kawamoto, H., Kanbe, S., & Sankai, Y. (2003). Power assist method for HAL-3 using
EMG-based feedback controller. In IEEE international conference on systems, manand cybernetics (pp. 1648–1653).

T. Matsubara et al. / Neural Networks 25 (2012) 191–199
Ko, J., & Fox, D. (2008). GP-Bayes filters: Bayesian filtering using Gaussian process
prediction and observation models. In IEEE/RSJ international conference on
intelligent robots and systems (pp. 3471–3476).
Lawrence, N. (2007). Learning for larger datasets with the Gaussian process latentvariable model. In Proceedings of the eleventh international workshop on artificial
intelligence and statistics (pp. 243–250).
Liu, X., & Goldsmith, A. (2004). Kalman filtering with partial observation losses. In
IEEE conferences on decision and control (pp. 4180–4186).
Li, Y., Wang, T., & Shum, H.-Y. (2002). Motion texture: a two-level statistical model
for character motion synthesis. ACM Transactions on Graphics, 21(3), 465–472.
Onishi, M., Luo, Z., Odashima, T., Hirano, S., Tahara, K., & Mukai, T. (2007).
Generation of human care behaviors by human-interactive robot RI-MAN. In
IEEE international conference on robotics and automation (pp. 3128–3129).
Ormoneit, D., Sidenbladh, H., Blank, M., & Hastie, T. (2001). Advances in neural
information processing systems: Vol. 13. Learning and tracking cyclic humanmotion (pp. 894–900).
Pavlovic, V., Rehg, J. M., & MacCormick, J. (2000). Advances in neural information
processing systems: Vol. 12. Learning switching linear models of human motion(pp. 981–987).
Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian processes for machine learning.
MIT Press.
Riley, M., Ude, A., Wada, K., & Atkeson, C.G. (2003). Enabling real-time full-body
imitation: a natural way of transferring human movement to humanoids. In
IEEE international conference robotics and automation (pp. 2368–2374).
Sato, M., & Ishii, S. (2000). On-line em algorithm for the normalized Gaussian
network. In Neural Computation, 12, 407–432.
Shapiro, A., Cao, Y., & Faloutsos, P. (2006). Style components. In Graphics interface(pp. 33–39).
Shumway, R. H., & Stoffer, D. S. (1982). An approach to time series smoothing andforecasting using the EM algorithm. Journal of Time Series Analysis, 3, 253–264.

199

Sidenbladh, H., Black, M. J., & Fleet, D. J. (2000). Stochastic tracking of 3D humanfigures using 2D image motion. In European conference computer vision. Vol. 2(pp. 702–718).
Taylor, G. W., & Hinton, G. E. (2009). Factored conditional restricted Bolzmann
machines for modeling motion style. In International conference on machinelearning (pp. 1025–1032).
Taylor, G. W., Hinton, G. E., & Roweis, S. T. (2006). Modeling human motion
using binary latent variables. In Proceedings of advances in neural information
processing systems (pp. 1345–1352).
Tenenbaum, J. B., & Freeman, W. T. (2000). Separating style and content with bilinear
models. Neural Computation, 12, 1247–1283.
Torresani, L., Hackney, P., & Bregler, C. (2006). Learning motion style synthesis fromperceptual observations. In Advances in Neural Information Processing Systems:
Vol. 19 (pp. 1393–1400).
Urtasun, R., Fleet, D.J., & Fua, P. (2006). 3D people tracking with Gaussian process
dynamical models. In IEEE Computer society conference on computer vision andpattern recognition (pp. 238–245).
Urtasun, R., Fleet, D.J., Hertzmann, A., & Fua, P. (2005). Priors for people tracking
from small training sets. In IEEE international conference on computer vision(pp. 403–410).
Urtasun, R., & Fua, P. (2004). 3D Human body tracking using deterministic
temporal motion models. In European Conference on Computer Vision. Vol. 3(pp. 92–106).
Wang, J. M., Fleet, D. J., & Hertzmann, A. (2007). Multifactor Gaussian process models
for style-content separation. In International conference on machine learning(pp. 975–982).
Wang, J. M., Fleet, D. J., & Hertzmann, A. (2008). Gaussian process dynamical models
for human motion. IEEE Transactions on Pattern Analysis and Machine Intelligence,30, 283–298.