Approximate Inference in Additive Factorial HMMs with
Application to Energy Disaggregation

J. Zico Kolter
MIT CSAIL
kolter@csail.mit.edu

Abstract
This paper considers additive factorial hidden Markov models, an extension to HMMs
where the state factors into multiple independent chains, and the output is an additive
function of all the hidden states. Althoughsuch models are very powerful, accurate inference is unfortunately difficult: exact inference is not computationally tractable, andexisting approximate inference techniques arehighly susceptible to local optima. In thispaper we propose an alternative inference
method for such models, which exploits theiradditive structure by 1) looking at the observed difference signal of the observation,2) incorporating a “robust” mixture component that can account for unmodeled observations, and 3) constraining the posterior toallow at most one hidden state to change at atime. Combining these elements we develop
a convex formulation of approximate inference that is computationally efficient, has
no issues of local optima, and which performs much better than existing approaches
in practice. The method is motivated by the
problem of energy disaggregation, the task oftaking a whole home electricity signal and decomposing it into its component appliances;applied to this task, our algorithm achieves
state-of-the-art performance, and is able toseparate many appliances almost perfectly
using just the total aggregate signal.

1

Introduction and Background

Factorial hidden Markov models (FHMMs) (Ghahramani and Jordan, 1997) are an extension of the basic
Appearing in Proceedings of the 15th International Conference on Artificial Intelligence and Statistics (AISTATS)
2012, La Palma, Canary Islands. Volume XX of JMLR:
W&CP XX. Copyright 2012 by the authors.

Tommi Jaakkola
MIT CSAIL
tommi@csail.mit.edu

hidden Markov model where several HMMs evolve independently in parallel, and the observed output issome joint function of all the hidden states. Due totheir ability to capture complex aggregate signals via
a compact structure, such models have found applications in areas as speech recognition (Virtanen, 2006),
audio separation (Roweis, 2001), and pitch tacking(Bach and Jordan, 2005). However, despite their utility, accurate inference is difficult: exact inference requires enumerating an exponential number of states,
and common approximate algorithms, (block) Gibbs
sampling (Kim et al., 2011) or structured mean field
methods (Ghahramani and Jordan, 1997), are highly
susceptible to local optima in the likelihood function.
In this paper we consider the special case of additive FHMMs, where each HMM emits an (unobserved)
real-valued output, and we observe the sum of these
outputs.1 The key algorithmic contribution of this paper is a method for exploiting the additive structure of
the FHMM to develop an approximate inference procedure that vastly outperforms existing approaches. Our
method considers models that represent both the total
aggregate output and the difference between successive
outputs, and also introduces a “robust” mixture component that can represent arbitrary additional signals(but with an ℓ1 -based regularity condition). By constraining the posterior to require that only one HMM
change state any any given time, we develop an efficient convex quadratic programming relaxation of theinference problem. The resulting method works extremely well in practice, is computationally efficient(scales roughly linearly in the number of HMMs), canhandle non-IID additive noise, and is free from local
optima in the optimization procedure.
We focus on the application of electrical energy disaggregation, also called non-intrusive load monitoring(Hart, 1992), the task of taking a whole-home energy
signal and breaking it down into its component appli1
The original FHMM model (Ghahramani and Jordan,
1997) was of this form, but since many other aggregation
functions have since been used and referred to as FHMMs,
we explicitly denote this case the additive FHMM.

Approximate Inference in Additive Factorial HMMs with Application to Energy Disaggregation

ances. Studies have shown that just presenting such a
breakdown to users, so that a home owner can see precisely how much energy is being used by which appliance, can automatically lead to energy-saving behavior (Darby, 2006; Neenan and Robinson, 2009). Given
the large percentage of energy that residential andcommercial buildings consume, enabling such energyconserving behavior is a key element to addressing energy sustainability problems. Although there has beenmuch work on energy disaggregation (see e.g. (Ziefman
and Roth, 2011) for a recent survey) it has largely focused simply upon classifying electrical events; early
work looked only at power signals along with finite
state machine models, (Hart, 1992), while later work
has incorporated transient and harmonic informationfrom very high-frequency sampling (Laughman et al.,
2003; Gupta et al., 2010; Berges et al., 2010). Recent work has begun to look at more complex inference
procedures for the actual disaggregation task, including FHMMs (with block Gibbs sampling) (Kim et al.,2011) and sparse coding methods (Kolter et al., 2010).
In comparison to past work, this paper makes several
contributions from an application standpoint. First,
unlike all past work we are aware of other than (Kimet al., 2011), we consider the unsupervised setting,where we do not have a priori knowledge of the appliances in the home; the inference procedure we develop, however, enables us to use much more complex
device models than this past work. Our method effectively joins the two prior threads of work in energy
disaggregation: we are able to capture general device
states and capture the information in transient powersignature, using a single unified model. Finally, ourinference method significantly outperforms alternative
inference methods for this task, and for several devices
is able to achieve essentially perfect separation.

2

The Additive Factorial Model

The basic FHMM consists of several independent
HMMs evolving in parallel, with the observation beinga joint function of all hidden states; a graphical model
representation is shown in Figure 1. Exact inference
typically is not tractable in such a model, so approximate inference procedures are needed, such as block
Gibbs sampling or structured mean field methods. In
this paper we are interested in the special case of the
FHMM where the output is an additive function of thedifferent hidden states as in (Ghahramani and Jordan,
1997). The factors in the additive FHMM are(i)x1 ∼ Mult(φ(i) )
(i)

(i)

(i)

xt |xt−1 ∼ Mult(P (i) )xt−1(1:N )
ȳt |xt
∼N

N
X
i=1

(i)

(1)
!

µ (i) , Σ
xt

...

xt−1

...

xt−1

xt

xt+1

..

..

..

ȳt−1

ȳt

ȳt+1

(1)

(2)

(1)

(1)

...

(2)

...

xt+1

xt

(2)

Figure 1: The FHMM model.
where N ∈ Z+ is the number of HMMs; mi ∈ Z+
is the number of states in the ith HMM; T ∈ Z+ is(i)the number of time steps; xt ∈ {1, . . . , mi } denotes
the state of the ith HMM at time t; ȳt ∈ Rn is the(i)observed aggregate output; µj ∈ Rn is the mean of
the ith HMM for state j; Σ ∈ Rn×n is the observation
variance; φ(i) ∈ [0, 1]mi is the initial state distribution
for ith HMM; and P (i) ∈ [0, 1]mi ×mi is the transition
matrix for ith HMM.
Although the additive model is a special case of thegeneral FHMM, exact inference is still not tractablein the model, and standard mean-field methods are
susceptible to local optima. To address these issues,
with an eye towards formulating a convex approximate
inference procedure, we propose two additions to themodel that guide the posterior.
2.1

The Difference FHMM ModelIn the additive FHMM above, it is natural to consider differences in the output signal: intuitively, inthe case that only one HMM changes state at a given
time, these differences correspond precisely to the statechange of this one HMM and can thus help us to infer
the hidden states. We will formalize this intuition in
the next section, but this has been a common heuristic in algorithms for energy disaggregation: looking for
sharp changes in the power signal and classifying theseas a single device change. We encode these difference
observations directly using a slight modification of theFHMM, which we refer to as the difference FHMM.
(i)
The model is shown in Figure 2; the the xt variables
and its factors are the same as for the FHMM, but theoutput at time t, denoted ∆ȳt , is given by
!
N 

X
(1:N )(1:N )∆ȳt |xt, xt−1 ∼ Nµx(i) − µx(i) , Σ
i=1

t

t−1

(2)
where we use the same parameters as for the additive
FHMM described earlier (we will hereafter use the abbreviation ∆µj,k = µj − µk ).
The additive FHMM and the difference FHMM defined here have complementary advantages and disadvantages. The additive FHMM captures well the total

J. Zico Kolter, Tommi Jaakkola

...

xt−1

...

xt−1

xt

xt+1

..

..

..

∆ȳt−1∆ȳt∆ȳt+1

(1)

(2)

(1)

(1)

...

(2)

...

xt+1

xt

(2)

directly represents the difference signal. We now introduce a signal ∆zt ∈ Rn , with prior
)
(
T
X
1
(5)
k∆zt k1
exp −λp(∆z1:T ) =Z(λ, T )t=1

and modify the likelihood of ∆ȳt to be(1:N )∆ȳt |xt

Figure 2: The difference factorial HMM model.

(1:N ),xt−1 , ∆zt ∼
N

N
X
i=1

aggregate output of the system, but does not directly
include information about the signal differences. In
contrast the difference FHMM directly encodes thesedifference observations, but fails to model the total aggregate output, so that errors in the cumulative sum of
this model can build up over time. Our final inference
procedure will thus ultimately use both models, butwill constrain the posteriors over the hidden states toagree between them.
2.2

Robust Mixture Components

The Gaussian likelihood term in the FHMM (both theadditive and difference versions) makes the model sensitive to outliers. This is an especially troublesome
problem for the energy disaggregation task, as there
may be new devices or rarely used devices for which
we have no model. We thus introduce a robust version of the additive mixture model, which includes a
generic component that can take on arbitrary values,
but which obeys certain regularity conditions. A particularly suitable choice of regularizer here is (1D) Total Variation (TV) Regularization (Rudin et al., 1992),
a penalty on the ℓ1 norm of the signal differences (orprobabilistically, a Laplace prior on the differences);
this prior encourages the generic mixture component
to take on piece-wise constant values, which intuitively
captures the typical nature of many devices.
To add this generic mixture component to the additive
FHMM model, we extend (1) by introducing a signal
zt ∈ Rn , with the prior
)
(
T
−1
X
1
kzt+1 − zt−1 k1 (3)exp −λp(z1:T ) =Z(λ, T )t=1
and modify the likelihood of the aggregate output to
be
N
X
(i)
(1:N )µ (i) + Σ1/2 zt , Σ)ȳt |xt, zt ∼ N (
(4)
i=1

xt

where we scale zt by Σ1/2 because we are using a single
penalty λ for all the dimensions of zt , and we want tomake its units commensurate with those of ȳt .
A similar robust version of the difference FHMM is
possible, and is even simpler in this case, as the model

1/2∆µx(i) ,x(i) + Σ
t

t−1

!

∆zt , Σ(6)

For intermediate values of λ (i.e.,a λ must be large
enough so that the model does not assign all the output
to zt but small enough so that z1:T is not a constantvalued), these additional mixture components are able
to explain elements of the signal that are not well modeled by any of the HMMs.

3

Approximate MAP Inference

For the FHMM model and the proposed extensions,exact inference remains a computationally expensive
procedure: we know of no way to perform exact inference short of enumerating all the states. Thus, inthis section, we focus on methods for approximate inference, and approximate MAP inference in particular.
By exploiting the additive structure of the FHMM, using both the additive and difference formulations, andby constraining the posterior, we develop a methodthat is computationally efficient, free of local optima,
robust to unknown components in the signal, andwhich in practice works much better than existing approaches.
3.1

Exact MAP Inference

We begin by considering optimization approaches toexact MAP inference. To define notation, we perform
optimization over the variables
o
n
(i)
(i)
(i)
Q = Q(xt ) ∈ Rmi , Q(xt−1 , xt ) ∈ Rmi ×mi . (7)These variables (in the exact MAP case) are in(i)(i)
dicators, so that Q(xt )j = 1 ⇔ xt = j and
(i)
(i)
(i)
(i)
Q(xt−1 , xt )j,k = 1 ⇔ xt−1 = j, xt = k. In thefollowing, all summations or constraints over t are
assumed to range from t = 1, . . . , T (or 2, . . . , T ,as appropriate), over i are assumed to range fromi = 1, . . . , N , and over j, k are assumed to range fromj, k = 1, . . . , mi . For the robust mixture components,
we also optimize over the variables zt , ∆zt ∈ Rn for
the additive and difference formulations respectively.
The first constraint on Q is that it must be locallyconsistent and must represent a valid distribution locally, (i.e., it must lie within the local marginal polytope (Wainwright and Jordan, 2008)) which for thisApproximate Inference in Additive Factorial HMMs with Application to Energy Disaggregation

problem takes the form


X
(i)


Q(xt )j = 1








j


X


(i)
(i)
(i)



Q(xt−1 , xt )j,k = Q(xt−1 )j 
L= Q: X
k


(i)
(i)
(i)



Q(xt−1 , xt )k,j = Q(xt )j 








k




(i)
(i)
(i)
Q(xt )j , Q(xt−1 , xt )k,j ≥ 0have little correspondence to the MAP solution. However, as we show in the next section, by constraining
the posterior in the difference FHMM formulation, we
can transform (10) to a mixed-integer linear program,which is much more amenable to convex relaxation.2(8)3.2

MAP inference can now be cast as optimizing the loglikelihood of the model, subject to the constraint thatQ ∈ L and that all the variables in Q take on binary
values (a constraint which we abbreviate Q ∈ {0, 1}).
Due to the Gaussian likelihood term, these problems
take the form of mixed-integer quadratic programs(MIQPs). For the additive FHMM model, this leads
to the optimization problem
minimize over {Q ∈ L ∩ {0, 1}, z1:T },
2

X (i)
1X
(i)
ȳt − Σ1/2 zt −µj Q(xt )j
2 t
i,j
Σ−1
X
(i)
(i)
(i)
+
Q(xt−1 , xt )j,k (− log Pk,j )(9)

As described intuitively above, a natural condition to
impose upon the posterior distribution over states isthat at most one HMM changes state at any given
time, a requirement we refer to as the one-at-a-time
condition. Indeed, if we consider each HMM in ourmodel to be a sufficiently discretized approximation
to a continuous time HMM, then this condition holdswith high probability.3 Formally,
Proposition 1. For each HMM i, let P (i) be the discretization of a continuous time Markov chain with incidence matrix Q(i) ∈ Rmi ×mi and discretization interval ∆t ∈ R+ . For total time tf ∈ R+ , this results in
an FHMM with tf /∆t time steps, and we let A denotethe event that there is some time where two HMMschange state simultaneously(i)t,i,j,k

+λ

X

(i)

(j)

A = {∃t, i 6= j : xt 6= xt+1 , xt

kzt − zt−1 k1 .

(j)

6= xt+1 }.

(11)

Then

t

p(A) ≤ O(∆t).

Similarly, exact MAP inference for the difference
FHMM model is given by the MIQP
minimize over {Q ∈ L ∩ {0, 1}, ∆z2:T },
2

X (i)
1X
(i)
(i)
∆ȳt − Σ1/2 ∆zt −∆µk,j Q(xt−1 , xt )j,k
2 t
i,j,k
Σ−1
X
X
(i)
(i)
(i)
+
Q(xt−1 , xt )j,k (− log Pk,j ) + λk∆zt k1 .
t,i,j,kThe One-at-a-time Condition and LP
Relaxations

t

(10)

Naturally, solving these problems exactly is difficult,as they involve non-convex integer constraints, andso some form of relaxation is necessary to obtain a
tractable formulation. A typical method for converting binary integer problems such as (9) and (10) into
tractable optimization problems is to drop the constraint that Q take on integer values, and solve theresulting, now convex, quadratic programs (this is
usually stated as relaxing that constraint that thevariables be in {0, 1} to the constraint that they bein [0, 1], but the later is already implied by the local marginal polytope in our case). Unfortunately,
the quadratic term in the optimization problems canmake this relaxation perform poorly in practice: this(i)term typically encourages the Q(xt ) variables to takeon non-integral values, and the resulting solution can(12)

The proof and definitions for a continuous time
Markov process are in Appendix A. Note the one-ata-time condition is not imposed directly by modifying
the FHMM model (for example, by adding a latentvariable that determines which of the HMMs change);
this would create an undesired dependence between all
the hidden states in subsequent time steps and complicate inference. Rather this condition is imposed asa constraint on the class of allowable posterior distributions, similar in spirit to the posterior regularization
framework (Ganchev et al., 2010), except that we arelooking only at the inference task rather than a joint
EM task, and the specific constraint we consider here
is of course quite different.
The set of distributions that obey the one-at-a-time
condition can be expressed formally as the constraint(1:N )(1:N )
that the off-diagonal elements of Q(xt−1 , xt
) sum
2

The same can be done in principle for the additivemodel, but would require us to introduce exponentiallymany marginal variables.
3
We could alternatively use a continuous time FHMM
model, where we assume that observations are generated(at least) whenever one of HMMs changes state; in this casethe one-at-a-time condition holds with probability one. We
use the discrete-time formulation in this paper for simplicity of the presentation.

J. Zico Kolter, Tommi Jaakkola

to at most one. That is,




X
(i)
(i)
O= Q:
Q(xt−1 , xt )j,k ≤ 1 .



(13)

i,j,k6=j

Imposing this constraint implies that at most one of
(i)
(i)
(i)
the ∆µk,j Q(xt−1 , xt )j,k terms in each norm in (10)is non-zero; further, if exactly one of these terms is
non-zero then ∆zt = 0, as we enforce the one-at-atime condition over the generic mixture component as
well, and so the quadratic term in (10) can be replaced
by a linear term. Alternatively, if none of the offdiagonal terms are non-zero then the quadratic term
is just k∆ȳt − Σ1/2 ∆zt k2Σ−1 , which (together with theλk∆zt k1 term) can be optimized analytically to givethe Huber loss function
D(y, λ) = min{ky − zk22 + λkzk1 }
z



n
X
λ2 λ2
1 2
y , max λ|yℓ | − ,
min
=
2 ℓ
2 2
ℓ=1

(14)

Thus, exact MAP inference for the difference FHMM
model (for the posterior constrained to the set O) canbe written as the mixed-integer linear program (seeAppendix B for the full derivation)
minimize over {Q ∈ L ∩ O ∩ {0, 1}},
2
1 X
(i)
(i)
(i)
∆ȳt − ∆µk,j −1 Q(xt−1 , xt )j,k2
Σt,i,j,k6=j


X
1X
(i)
(i)
D(Σ−1/2 ∆ȳt , λ) 1 −Q(xt−1 , xt )j,j 
+
2 t
i,j
X
(i)
(i)
(i)
+
Q(xt , xt+1 )j,k (− log Pk,j )t,i,j,k(15)

Dropping the integral constraint in (15) leads to a convex LP. While solutions to the LP need not be integer
valued, in practice they often are: the objective termand simplex constraints in (15) alone result in separable optimization problems that would always haveinteger solution, and the only constraints that may
cause non-integral solutions are the consistency terms
in the local marginal polytope, which typically leads
to a very small number of non-integral assignments.
3.3

Joint Approximate Inference

At this point we have two potential models for performing inference which capture different elements ofthe signal: the total aggregate output and the difference signal. To exploit the benefits of each model
we combine these optimization problems (after relaxing the integer constraint) into a single joint problem;

Input: ȳ1:T ∈ Rn , aggregate output signal;(1:N )µ1:mi ∈ Rn , state means for N HMMs; Σ1 , Σ2 ∈Rn×n , λ1 , λ2 ∈ R+ , covariance and regularization
parameters
minimize over {Q ∈ L ∩ O, z1:T }
2

X (i)
1X
1/2
(i)
ȳt − Σ1 zt −µj Q(xt )j
2 t
i,j
+

1
2

X

(i)

∆ȳt − ∆µk,j

t,i,j,k6=j

2
Σ−1
2

(i)

Σ−1
1
(i)

Q(xt−1 , xt )j,k



X
1X
−1/2
(i)
(i)
+
D(Σ2 ∆ȳt , λ2 )1 −Q(xt−1 , xt )j,j
2 t
i,j
X
(i)
(i)
(i)
+
Q(xt−1 , xt )j,k (− log Pk,j )t,i,j,k+ λ1

Xkzt − zt−1 k1t(1:N )

Output: ŷ1:T , predicted individual HMM output
X (i)
(i)
(i)
ŷt =
µj Q(xt )jj

Figure 3: The AFAMAP algorithm.
equivalently, we are performing inference using bothmodels, but constrain the posterior over the hidden
variables in both cases to agree (note that we don’t
require the zt variables to agree, but rather they act
to make each model robust to outliers). This leads
to the final optimization problem shown in Figure 3,
and we refer to the resulting algorithm as AFAMAP(Additive Factorial Approximate MAP).
The AFAMAP optimization problem is a convex
quadratic program, so could in principle be solved using off-the-shelf libraries. However, the total variation term in particular is poorly handled by generalsolvers. Fixing Q, however, we are left with a simple proximal TV regularization problem, which is efficiently solved using recent methods (Barbero and Sra,2011). Thus, we employ alternating minimization, using a generic QP solver (we use the CPLEX solver in
our implementation) to solve the AFAMAP problem
for fixed z1:T , then using the custom solver referenced
above to find z1:T given a fixed Q. A major advantage to the AFAMAP algorithm is that the objective
is jointly convex Q and z1:T and so this procedure willfind the global optimum; if we try to incorporate arobust mixture component using standard MAP inference, some variational approximation would be necessary. Also crucial to obtaining good performance is toproperly exploit sparsity: if a transition matrix P (i) is

Normalized Disaggregation Error

Approximate Inference in Additive Factorial HMMs with Application to Energy Disaggregation

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1

Normalized Disaggregation Error

0

2

4

6
Number of HMMs

8

10

8

10

Exact MAP
AFAMAP
Structured Mean Field
AFAMAP, no diff
AFAMAP, no add
Prior

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

Normalized Disaggregation Error

Exact MAP
AFAMAP
Structured Mean Field
AFAMAP, no diff
AFAMAP, no add
Prior

2

4

6
Number of States

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

2

4
6
8
Dimensionality of Output

2t,iExact MAP
AFAMAP
Structured Mean Field
AFAMAP, no diff
AFAMAP, no add
Prior

0.8

set.4 For state spaces that are small enough we run exact MAP inference. We evaluate the above-mentioned
structured variational mean field (SMF) algorithm(Ghahramani and Jordan, 1997), which approximates
the FHMM posterior by N decoupled HMMs, and iteratively improves the log-likelihood. We also evaluate the AFAMAP algorithm with one of the two components removed (i.e., with Σ1 = ∞ or Σ2 = ∞).
While the objective of SMF is slightly different, as it iscomputing a full distribution over hidden states rather
than just a MAP estimate, in virtually all the cases weconsider (both here and in the next section) the posterior is sufficiently peaked that the outputs are usually
nearly integer valued. Finally, as a baseline we compare to a method that simply sets the expectation of
all states to the stationary distribution of the Markov
chain. To evaluate the methods, we use the normalized disaggregation error, which directly measures how
well the models recovered the individual HMM out(i)(i)
puts: given true output yt and predicted output ŷt ,this is defined as
v
 

u
u X
X
2
2
u
(i)
(i)
(i)
/

t
yt − ŷt
yt
(16)

10

Figure 4: Performance varying the number of HMMs
(top), the number of states (middle), and the dimensionality of the output (bottom).
sparse, then we only need to represent those values of(i)(i)
Q(i) (xt−1 , xt ) corresponding to non-zero entries.
Given a (possibly non-integral) solution of the optimization problem, we could potentially refine it using
branch and bound approaches, randomized rounding,or run exact MAP inference only over the non-integral
states. In practice, however, since the LP objective
encourages integer-valued solutions, simply predicting
the mean output works well, and we use this approachin AFAMAP.

Experiments on Synthetic Data

To measure the effectiveness of our inference algorithmwith respect to different elements of the problem setting, we evaluate the AFAMAP algorithm versus several potential competitors on a simple synthetic data

2

Figure 4 shows the performance of the different algorithms varying the number of HMMs, states per
HMMs, and the dimensionality of the output. In
all cases AFAMAP outperforms all alternative approaches (except for exact MAP inference when applicable). The failure of the SMF approach is particularly
evident here: because the algorithm assigns states to
HMMs one at a time, it often will explain some portion of the aggregate signal using the incorrect state;even using annealing-based approaches (inflating thevariance in early iterations, which we do for all SMF
experiments), it is often difficult to get out of this local
optimum. The AFAMAP algorithm, in contrast, is aconvex method and thus has no problems of local optima; furthermore, training using both the difference
and additive models outperforms either in isolation.
Our second set of experiments evaluates the robust
mixture component. Here we add a random walk
4
Experimental details: To create the data, we generate N “cyclic” HMMs, each with m states, and an ndimensional output; the initial state distribution is uniform over all states, the transition matrix is generated by(i)(i)
Pj,j ∝ 30, Pmod(j+1,k+1)+1,j ∝ U [1, 2] and the mean of
(i)

4

t,i

each HMM is (µj )ℓ ∼ U [0, 2]. For each HMM we sample
T = 500 time steps, and observe the sum of all the individual outputs plus Gaussian noise sampled from N (0, 0.01I).
For all experiments we use N = 4, m = 4, and n = 4,
except if one of these parameters is being varied in the experiment. Regularization parameters were set to be on theorder of the parameters of the true model. MATLAB codefor the experiments are included with the paper.

Normalized Disaggregation Error

J. Zico Kolter, Tommi Jaakkola

microwave

AFAMAP
Exact MAP
AFAMAP, Robust
Exact Map, Robust
Prior

0.8
0.7
0.6

1500

1000
500

0

0.4

0

5

10 15

0

500
0

bath gfi

0.3
0.2

1500

0.1

1000

100

200

0

0

electronics

20

40

60

kitchen outlet
1000

1000

500

500

500

0

0.02
0.04
0.06
0.08
Std Dev. of Random Walk Noise

0.1

0

0

100

200

signal to the output, ȳt = ȳt + wt where wt ∼N (wt−1 ; σ 2 I), and varied σ from 0 to 0.1. Figure
5 shows the performance of “exact” MAP and the
AFAMAP algorithm, with and without the robust
component. This experiment highlights the benefit of
the convex formulation: for the exact MAP formulation, we need some way of including the TV mixture(1:N )
component, since joint MAP inference over x1:T andz1:T is not tractable; in this example, we used simple
alternating maximization which leads to similar issues
as with SMF: initial iterations “lock” the HMM into
bad state guesses, and the resulting optimization optimization over the robust mixture cannot correct this.
Alternatively, the convex objective of AFAMAP is able
to estimate the noise component.

Electrical Energy Disaggregation

As mentioned above, the primary application we areconcerned with in this paper is the task of electrical energy disaggregation, separating a whole-home energy
signal into its component appliances. Although the focus of this paper is the inference procedure, briefly, thesetup is as follows. Using a high-resolution A/D device, we logged whole-home aggregate power consumption for two weeks. We also logged power consumption at the circuit level to obtain ground truth labels(this was used only for validation, and not for training the models). To reduce the signal size, we useda greedy variant of total variation regularization toapproximate the power signal as piecewise constant.5
Using the total power signal, we extract all snippets ofdata where consumption increases over some threshold then eventually returns to its original level (i.e.,
where some device comes on and turns off); some ofthese are indecipherable due to the aggregation, but
especially for short device durations there are occasional snippets of individual devices. We model allthese snippets as “empirical” HMMs (means equal to
All data is available at http://redd.csail.mit.edu.

0furnace

0

20

40

0

0

kitchen outlet

10

20

washer dryer

400
1000

2000
200

500

5

1500

1000

500

Figure 5: Performance of algorithms, with and withoutrobust mixture component, with random walk noise ofdifferent levels added to the output.

5bath gfi

1000

0.5

0

bath gfi

00

200 400 600

01000
0 20 40 60

0

0

50

100

Figure 6: Extracted signatures for nine appliancetypes, labeled by device circuit, with time in seconds
on the x-axis and power on the y-axis. Each stateof the device HMM corresponds to a constant-power
segment in the plot, with length equal to that state’s
expected duration.
the observed power output, and transitions probabilities set based upon the amount of time spent at eachpower level) and looked at all pairwise probabilities
between them (the probability of one snippet generating another); using the k-nearest-neighbor graph induced by these probabilities, we ran spectral clustering
to group devices together. This resulted in selecting
nine “prototypical” motifs, shown in Figure 6, that
occurred frequently in the data, and which correspondvery clearly to different devices.
Treating these nine extracted motifs as the HMMs in afactorial model, we ran the AFAMAP to separate theappliances in the aggregate power signal over the entire
two week period. For all methods, regularization parameters (λ and Σ) were fit using one day of the data.
Table 1 shows the performance of the methods for eachdevice. Because we do not know the true contributionof the device (just the circuit contribution), we report
precision and recall metrics at the circuit level: recall measures what portion of a given circuit’s energy
is correctly classified, while precision measures, of theenergy assigned to a circuit, how much truly belonged
to that circuit. Because there are some portions of theenergy that are correctly left unassigned (because wehave no model of that device), we would not expect toachieve 100% recall, and we place higher emphasis on
the precision metric.
For all the circuits in the home, AFAMAP vastly outperforms the SMF method on this problem. In particular, for five of the nine devices considered, AFAMAPApproximate Inference in Additive Factorial HMMs with Application to Energy Disaggregation

5000

Power (Watts)

4000

Table 1: Per-circuit performance for AFAMAP and
SMF for two weeks of data. Performance is reported
as precision/recall, and bold entries denote staticallysignificant better performance on both metrics. Morethan one device above can belong to a single circuit.

kitchen outlets
microwave
electronics
washer dryer

3000
2000
1000
0

17:40

17:50

18:00
Time

18:10

18:00
Time

18:10

18:00
Time

18:10

5000

Power (Watts)

4000

unassigned
kitchen outlets
microwave
washer dryer

3000

Circuit
1 Microwave
2 Bath GFI
3 Electronics
4 Kitch. Out. 1
5 Furnace
6 Kitch. Out. 2
7 Wash/Dryer
Total

AFAMAP
97.5% / 66.1%
82.7% / 70.8%
41.6% / 0.8%
37.5% / 12.9%
91.7% / 70.8%
45.2% / 16.0%
98.8% / 73.6%
87.2% / 60.3%

SMF
96.8% / 4.1%
50.1% / 9.1%
41.4% / 0.3%
10.2% / 47.8%
12.6% / 15.3%
13.3% / 24.8%
89.3% / 76.7%
35.5% / 45.1%

2000
1000
0

Power (Watts)

5000
4000
3000
2000

17:40

17:50

unassigned
bath gfi
microwave
furnace
kitchen outlets
kitchen outlets
washer dryer

robust TV mixture component to capture the remaining power, for which it does not have a model. In contrast, SMF assigns the power to the washer/dryer (almost) correctly, but completely fails to properly assign
the other devices. This is indicative of the method’s
tendency to “lock in” to variable assignments; indeed,
we tuned SMF substantially to get even this result: optimizing the HMMs in a random order, for example,
typically misses the washer/dryer entirely.

1000
0

17:40

17:50

Figure 7: True (top), AFAMAP predicted (middle),
and SMF predicted (bottom) breakdown by device fora short portion of power data.
has a precision in its assignments of more than 90%;
this starts to reach the order of the accuracy in sensing between the two data sources, and a visual inspection of the signals confirms that AFAMAP is assigning these devices essentially perfectly. AFAMAP does
make errors in the remaining four devices, which uponinspection is due to the fact that appliances on multiple circuits have signatures very similar to these; thefive devices identified correctly have distinct patterns
that AFAMAP can identify even when they occur inaggregated form. In contrast, the SMF algorithm essentially fails completely on this task: either the algorithm predicts a great deal of energy incorrectly, or it
assigns very little energy to a circuit.
The difference between the AFAMAP and SMF algorithms is equally striking if we look visually at thepredicted separations of the two methods. Figure 7
looks at a short period of time and shows the true
appliance breakdown along with the predictions of
AFAMAP and SMF. AFAMAP correctly identifies thewasher/dryer, microwave, and kitchen outlets during
this time (with minor errors), and (correctly) uses the6

Conclusion

Although the additive factorial HMM model has large
representation power, its applicability has likely been
greatly reduced by the difficulty of inference in models with a substantial number of HMMs. By exploiting the additive structure of the problem and by constraining the set of allowed posteriors, the AFAMAP
algorithm can accurately perform inference on suchtasks, opening the door to many more potential applications. While this paper focused specifically on
the inference problem and the application to energy
disaggregation, it would be natural to investigate the
combination of learning and inference using this procedure, as well as look at using the method as an initialization for other inference procedures. From anapplication standpoint, many extensions look promising as well. The unsupervised learning procedure we
described briefly is currently limited to finding devices
with relatively short time scales; but by applying thisapproximate inference and iteratively looking at the“unassigned” portions of energy, it may be possible tosuccessively build models for more and more devices
in the home. Combined with joint inference and “hardEM” learning procedures (EM but using MAP inference) these methods have the potential to achieve theeventual goal of disaggregating virtually all the appliances in a home without supervision.

J. Zico Kolter, Tommi Jaakkola

References
Bach, F. R. and Jordan, M. I. (2005). Discriminative
training of hidden markov models for multiple pitch
tracking. In Proceedings of the International Conference on Accoustics, Speech, and Signal Processing.
Barbero, A. and Sra, S. (2011). Fast newton-type
methods for total variation regularization. In Proceedings of the International Conference on Machine
Learning.
Berges, M., Goldman, E., Matthews, H. S., and Soibelman, L. (2010). Enhancing electricity audits in residential buildings with non-intrusive load monitoring. Journal of Industrial Ecology: Special Issue
on Environmental Applications of Information andCommunications Technology, 14(5):844–858.
Darby, S. (2006). The effectiveness of feedback on energy consumption. Technical report, Environmental
Change Institute, University of Oxford.
Ganchev, K., Graça, J., Gillenwater, J., and Taskar,
B. (2010). Posterior regularization for structured
latent variable models. Journal of Machine Learning
Research.
Ghahramani, Z. and Jordan, M. I. (1997). Factorial
hidden markov models. Machine Learning, 27:245–
273.
Gupta, S., Reynolds, S., and Patel, S. N. (2010). ElectriSense: Single-point sensing using EMI for electrical event detection and classification in the home. In
Proceedings of the Conference on Ubiquitous Computing.
Hart, G. (1992). Nonintrusive appliance load monitoring. Proceedings of the IEEE, 80(12).
Kim, H., Marwah, M., Arlitt, M., Lyon, G., and Han,
J. (2011). Unsupervised disaggregation of low frequency power measurements. In Proceedings of the
SIAM Conference on Data Mining.
Kolter, J. Z., Batra, S., and Ng, A. Y. (2010). Energy
disaggregation via discriminative sparse coding. In
Neural Information Processing Systems.
Laughman, C., Leeb, S., and Lee (2003). Advanced
non-intrusive monitoring of electric loads. IEEE
Power and Energy.
Neenan, B. and Robinson, J. (2009). Residential electricity use feedback: A research synthesis and economic framework. Technical report, Electric Power
Research Institute.
Roweis, S. T. (2001). One microphone source separation. In Neural Information Processing Systems.
Rudin, L. I., Osher, S., and Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms.
Phyica D, 60:259–268.

Virtanen, T. (2006). Speech recognition using factorial
hidden markov models for separation in the feature
space. In Procedings of Interspeech.
Wainwright, M. J. and Jordan, M. I. (2008). Graphicalmodels, exponential families, and variational inference. Foundations and Trends in Machine Learning,
1(1–2):1–305.
Ziefman, M. and Roth, K. (2011). Nonintrusive appliance load monitoring: Review and outlook. IEEE
Transactions on Consumer Electronics, 57(1):76–84.

A

Proof of Proposition 1

A continuous time Markov chain with m states is defined via an incidence matrix Q ∈ Rm×m . Assuming
the system is is state j at real-valued time t ∈ R+ ,we sample a duration from an exponential distribution with rate −Qjj . We then transition to statex(t + dt) = k with probability
P −Qk,j /Qj,j (which implies that Q must have
k6=j Qk,j = −Qj,j ). Thatisdt|x(t) = j ∼ Exp(−Qj,j )x(t + dt) ∼ Mult(q), qk =

(

−Qk,j
Qk,k

0

k 6= jk=j(17)

For some time interval ∆t, the discretized transition
matrix, defined asP (∆t)k,j = p(x(t + ∆t) = k|x(t) = j),(18)

is given byP (∆t) = exp(∆tQ)(19)

where exp(·) here is the matrix exponential defined asexp(X) =

∞
X
1

i!
i=0

Xi

(20)

We now prove Proposition 1 using these definitions.
Proof. We first note that for r = − minj Qj,j (thelargest rate parameter for the exponential distributions)min P (∆t)j,j ≥ 1 − r∆t.
(21)jThis follows immediately from the matrix inequality
exp(X)  I + X which in turn follows from the eigenvalue representation of matrix exponentiation
exp(X) = U exp(Λ)U −1(22)where X = U ΛU −1 is the eigenvalue decompositionof X. Since Λ is diagonal, the matrix exponential issimply the elementwise exponentiation of each entry
on the diagonal, and the matrix inequality follows by

Approximate Inference in Additive Factorial HMMs with Application to Energy Disaggregation

applying the basic inequality ex ≥ 1 + x for each of the
diagonal entries.
Now, for a single time step let A0 be the event that
only zero or one of the N HMMs change state. Us(i)ing the bound above with r = p − mini,j Qj,j definedover all the HMMs with ∆t sufficiently small such thatr∆t < 1,p(A0 ) ≥ (1 − r∆t)N + N r∆t(1 − r∆t)N −1≥ (1 − r∆t)N + N r∆t(1 − r∆t)N(23)= (1 + N r∆t)(1 − r∆t)N .
Then probability of the one-at-a-time condition holding for all tf /∆t time steps is just P (¬A) ≥P (A0 )tf /∆t andby considering the mixed-integer quadratic program(10) along with the one-at-a-time constraint (13). As
mentioned in the text, enforcing the one-at-a-time constraint implies that at most one term in the summationX(25)

log(1 + N r∆t) + N log(1 − r∆t)
1
1
= N r∆t − N 2 r2 ∆t2 − N r∆t − N r2 ∆t2 + O(∆t3 )22
1= − (N + N 2 )r2 ∆t2 + O(∆t3 ).
2(26)
Combining this with (24) gives
1
log P (¬A) ≥ − tf (N + N 2 )r2 ∆t + O(∆t2 ). (27)2Finally, using the inequality ex ≥ 1 + x, we have
P (A) = 1 − P (¬A)
= 1 − exp(log P (¬A))(28)≤ − log P (¬A)
1
≤ tf (N + N 2 )r2 ∆t + O(∆t2 ) = O(∆t)
2

(i)

(i)

(29)

i,j,k

is non-zero.
First suppose there is one HMM i at time t with an(i)(i)
off-diagonal entry Q(xt−1 , xt )j,k , j 6= k. Then theerror in approximating ∆ȳt will be
(i)

(i)

(i)

∆ȳt − ∆µk,j Q(xt−1 , xt )j,ktf(log(1 + N r∆t) + N log(1 − r∆t)) .
log P (¬A) ≥∆t(24)
Using the Taylor approximation1log(1 + x) = x − x2 + O(x3 )2
we have that(i)∆µk,j Q(xt−1 , xt )j,k(i)∆ȳt − ∆µk,j

=

(i)

2
Σ

2
Σ−1

(30)

(i)
(i)
Q(xt−1 , xt )j,k−1(i)since Q(xt−1 , xt )j,k = 1 by the integrality constraint,(i′ )and since ∆µj ′ ,j ′ = 0 for all other (diagonal) entries
(i′ )

(i′ )

Q(xt−1 , xt )j ′ ,j ′ = 1. Furthermore, by the one-at-a(i)(i)
time condition, the ∆zt = 0 when Q(xt−1 , xt )j,k = 1.
Thus, the quadratic error term
2

X
1
(i)
(i)
(i)
∆ȳt − Σ−1/2 ∆zt −∆µk,j Q(xt−1 , xt )j,k
2
i,j,k

Σ−1

(31)

is equivalent in this case to the linear error term
1 X
2i,j,k6=j(i)∆ȳt − ∆µk,j

2
Σ−1

(i)

(i)

Q(xt−1 , xt )j,k

(32)

as desired, and the ℓ1 term k∆zt k1 = 0, so does notaffect the optimization problem.
Alternatively, suppose that no HM M changes state attime t. Then the error in approximating ∆ȳt is simply

While we focus on the dependence on ∆t in our presentation of the theorem, it is worth noting that the dependence is also quadratic in the number of HMMs, N .
Thus, for large numbers of concurrent HMMs (where
the probability of changing state is not overly small),
we may still need a fairly high sampling frequency.

B

Derivation of the MILP for
Difference FHMMs

Here we present a detailed description of how we arrive at the mixed-integer linear programming formulation to MAP inference in the difference FHMM (15),11k∆ȳt − Σ1/2 ∆zt k2Σ−1 = kΣ−1/2 ∆ȳt − ∆zt k22 (33)22This term plus the ℓ1 penalty on ∆zt can be optimized analytically, and has the well-known “softthresholding” solution,1
arg min ky − zk22 + λkzk1 = sign(y) max{|y| − λ, 0}z 2(34)
where all the operations in the right hand side areapplied elementwise to the entries of y. This solution

J. Zico Kolter, Tommi Jaakkola

attains objective value
n
X
1

2
i=1

(yi − sign(yi ) max{|yi | − λ, 0})2 + λ max{yi − λ, 0}

n
X

1
1
min{ yi2 , λ2 } + λ max{|yi | − λ, 0}
2
2
i=1



n
X
1 2
λ2 λ2
min
=
y , max λ|yi | − ,
2 i
2 2
i=1
=

= D(y, λ)(35)
Thus, the error (31) is equivalent in this case to theerror term


X
(i)
(i)
−1/2
D(Σ
∆ȳt ) 1 −Q(xt−1 , xt )j,k  (36)i,j,k6=j

P
(i)
(i)
since
i,j,k6=j Q(xt−1 , xt )j,k = 0 when no HMMschange state. Combining (32) and (36) give the MILP
formulation (15) as desired, and have eliminated the∆zt variable from the optimization problem.